\documentclass[preprint,aps,pra,onecolumn,superscriptaddress]{revtex4-1} %reprint
%\tightenlines

%\draft
\usepackage{etex}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{bbm}
\usepackage{listings}
% % \textwidth 16cm \textheight 23.5cm
% \renewcommand{\baselinestretch}{1.2}
\usepackage{graphicx}
\usepackage{graphics}
\usepackage{epsfig}
\usepackage{color}
\usepackage[dvipsnames]{xcolor}
\usepackage{multirow}
\usepackage[colorlinks]{hyperref}
\usepackage{fancyhdr}
\usepackage{calc}
\usepackage{natbib} %[numbers]
\usepackage{bibentry}
\usepackage{bbm}
\usepackage{siunitx}

% todo list and commands
%\usepackage{todonotes}
%% to avoid the conflict with amths package % not working
%\makeatletter
%\providecommand\@dotsep{5}
%\makeatother
%\listoftodos\relax
%\usepackage{makeidx}
%\allowdisplaybreaks
%% for eps transfering to pdf.
%\usepackage[update,prepend]{epstopdf}
%\usepackage{ifpdf}
%
%\ifpdf
%   \usepackage{graphicx}
%   \usepackage{epstopdf}
%   \epstopdfsetup{suffix=}
%   \DeclareGraphicsRule{.eps}{pdf}{.pdf}{`epstopdf #1}
%   \pdfcompresslevel=9
%\else
%   \usepackage{graphicx}
%\fi
% subfig
\usepackage{mwe}
\usepackage[caption=false]{subfig}
% to fix a figure's position using [H] option of thec figure.
\usepackage{float}
% to use \lesssim and other math symbols
\usepackage{amssymb}
% To include tikz plots.
\usepackage{tikz,pgfplots}
% recommended as of Pgfplots 1.3:
%\pgfplotsset{compat=newest}
\pgfplotsset{ % Here we specify options for all figures in the document
  compat=1.8, % Which version of pgfplots do we want to use?
  legend style =
  {font=\small\sffamily},
  label style = {font=\small\sffamily},
  xticklabel style={/pgf/number format/.cd, 1000 sep={}}
  }
%\usepackage[abs]{overpic}
% Precompile PDF figures from Tikz. Code is adapted from http://neishabouri.net/tips/2014/01/27/tips-and-tricks-for-creating-figures-in-tikz.html and the tex file has to be compiled with pdflatex -shell-escape file.tex format.
\usetikzlibrary{external} \tikzset{external/system call={lualatex
    \tikzexternalcheckshellescape -halt-on-error -interaction=batchmode -jobname "\image" "\texsource"}}
\tikzexternalize[prefix=fig/autofig/]


% self-defined short-cuts and commands
%\input{Mydef.tex}
\DeclareMathOperator{\tr}{tr}
\newcommand{\dt}[1]{\frac{{\mathrm d} {#1}}{{\mathrm d}t}}
\def\br{\mathbf{r}}
\def\bra#1{\langle{#1}\rvert}%{\mathinner{\langle{#1}\rvert}}
\def\ket#1{\lvert{#1}\rangle}%{\mathinner{\lvert{#1}\rangle}}
\def\Braket#1#2{\mathinner{\langle{#1}\! \mid\! {#2} \rangle}}
%========================================================================================
\newcommand{\erf}[1]{Eq.~(\ref{#1})}
\newcommand{\frf}[1]{Fig.~\ref{#1}}
\newcommand{\srf}[1]{Sec.~\ref{#1}}
\newcommand{\nn}{\nonumber}
\newcommand{\mbf}[1]{\mathbf{#1}}
%========================================================================================
% General quantum mechanics macros
%========================================================================================
\newcommand{\op}[2]{\ket{#1}\bra{#2}}
\newcommand{\expt}[1]{\langle{#1}\rangle}
\newcommand{\dg}{^\dagger}
\newcommand{\smallfrac}[2]{\mbox{$\frac{#1}{#2}$}}
\newcommand{\Tr}{\mbox{Tr}}
%========================================================================================
\newcommand{\expect}[1]{\big\langle #1 \big\rangle}
\newcommand{\eff}{\text{eff}}



% Redefine the tensor command.
%\renewcommand{\tensor}[1]{\boldsymbol{#1}}


%==== Ben's new macros ======
%\newcommand{\srf}[1]{Sec. \ref{#1}}
\newcommand{\half}{\smallfrac{1}{2}}

%==== subscripts ======
\newcommand{\oneD}{{\rm 1D}}
\newcommand{\vac}{{\rm vac}}
\newcommand{\cav}{{\rm cav}}
\newcommand{\inp}{{\rm in}}
\newcommand{\out}{{\rm out}}
\newcommand{\inter}{{\rm int}}
\newcommand{\scs}{{\rm SCS}}
\newcommand{\fwd}{+}
\newcommand{\bwd}{-}
\newcommand{\trans}{+}
\newcommand{\refl}{-}

 %==== operators/moments ======
\newcommand{\der}[1]{\frac{d {#1}}{dt}}
\newcommand{\unittens}{\tensor{\mathbf{I}}}
\newcommand{\poltens}{\hat{\tensor{\boldsymbol{\alpha}}}}
\newcommand{\varz}{\Delta J_3^2}
\newcommand{\jx}{\hat{J}_1}
\newcommand{\jz}{\hat{J}_3}
\newcommand{\shotnoise}{\Delta \mathcal{M}^2 |_{\rm SN}}
\newcommand{\projnoise}{\Delta \mathcal{M}^2_{\rm PN}}
\newcommand{\polcomp}{\hat{K}} % p,p' component of the tensor polarizability
\newcommand{\fo}{\hat{\mathbf{f}}}
\newcommand{\fx}{\hat{f}_x}
\newcommand{\fy}{\hat{f}_y}
\newcommand{\fz}{\hat{f}_z}
\newcommand{\Fx}{\hat{F}_x}
\newcommand{\Fy}{\hat{F}_y}
\newcommand{\Fz}{\hat{F}_z}
\newcommand{\rhoo}{\hat{\rho}}

%==== Microscopic moments for the qubit/qutrit subspace ====
\newcommand{\sigmauu}{\hat{\sigma}_{\uparrow\uparrow}}
\newcommand{\sigmaud}{\hat{\sigma}_{\uparrow\downarrow}}
\newcommand{\sigmaut}{\hat{\sigma}_{\uparrow \mathrm{T}}}
\newcommand{\sigmadu}{\hat{\sigma}_{\downarrow\uparrow}}
\newcommand{\sigmadd}{\hat{\sigma}_{\downarrow\downarrow}}
\newcommand{\sigmadt}{\hat{\sigma}_{\downarrow \mathrm{T}}}
\newcommand{\sigmatu}{\hat{\sigma}_{\mathrm{T}\uparrow}}
\newcommand{\sigmatd}{\hat{\sigma}_{\mathrm{T}\downarrow}}
\newcommand{\sigmatt}{\hat{\sigma}_{\mathrm{T}\mathrm{T}}}
\newcommand{\sigmaab}{\hat{\sigma}_{ab}}
\newcommand{\sigmaba}{\hat{\sigma}_{ba}}
\newcommand{\sigmadc}{\hat{\sigma}_{dc}}
\newcommand{\sigmacd}{\hat{\sigma}_{cd}}

\newcommand{\Dsigmauu}{\Delta\sigma_{\uparrow\uparrow}}
\newcommand{\Dsigmaud}{\Delta\sigma_{\uparrow\downarrow}}
\newcommand{\Dsigmaut}{\Delta\sigma_{\uparrow \mathrm{T}}}
\newcommand{\Dsigmadu}{\Delta\sigma_{\downarrow\uparrow}}
\newcommand{\Dsigmadd}{\Delta\sigma_{\downarrow\downarrow}}
\newcommand{\Dsigmadt}{\Delta\sigma_{\downarrow \mathrm{T}}}
\newcommand{\Dsigmatu}{\Delta\sigma_{\mathrm{T}\uparrow}}
\newcommand{\Dsigmatd}{\Delta\sigma_{\mathrm{T}\downarrow}}
\newcommand{\Dsigmatt}{\Delta\sigma_{\mathrm{T}\mathrm{T}}}
\newcommand{\Dsigmaab}{\Delta\sigma_{ab}}
\newcommand{\Dsigmaba}{\Delta\sigma_{ba}}
\newcommand{\Dsigmadc}{\Delta\sigma_{dc}}
\newcommand{\Dsigmacd}{\Delta\sigma_{cd}}

%==== physical parameters ======
\newcommand{\Eamp}{\mathcal{F}_0^{(+)}}
\newcommand{\charpol}{\alpha_0(\Delta_{f\!f'})}
\newcommand{\charpolq}{\alpha_0(\Delta_{f\!f'}^q)}
\newcommand{\qaxis}{\mathbf{e}_{\tilde{z}}}
\newcommand{\qangle}{\varphi}
\newcommand{\magic}[1]{\tilde{\omega}_{#1}}
\newcommand{\chiN}{\chi_{N}}
\newcommand{\NA}{N_C}
\newcommand{\chieff}{\chi_{\raisebox{-.1pt}{\tiny $J_3$}}}

%==== scattering and optical pumping rates ====%
\newcommand{\gammauu}{\gamma_{\uparrow \rightarrow \uparrow}}
\newcommand{\gammadd}{\gamma_{\downarrow \rightarrow \downarrow}}
\newcommand{\gammaud}{\gamma_{\uparrow \rightarrow \downarrow}}
\newcommand{\gammadu}{\gamma_{\downarrow \rightarrow \uparrow}}
\newcommand{\gammau}{\gamma_{\uparrow}}
\newcommand{\gammad}{\gamma_{\downarrow}}

%==== effective areas ======
\newcommand{\Ain}{A_{\rm in}}
\newcommand{\Abir}{A_N}
\newcommand{\AF}{A_{\rm Far}} % for the Faraday protocol.
\newcommand{\Ai}{A_0} % for the input light.
\newcommand{\Aint}{A_{\rm int}} % for the interaction area.

%==== eigenfunctions ======
\newcommand{\eigenf}{\mbf{f}_\eta}
\newcommand{\eigenfp}{\mbf{f}_{\eta'}}
\newcommand{\eigeng}{\mbf{g}_\eta}
\newcommand{\eigengp}{\mbf{g}_{\eta'}}

%==== field operators ======
\newcommand{\awg}{\hat{a}_{b,p}(\omega)}
\newcommand{\awr}{\hat{a}_{m,p}(\omega,\beta)}

%==== Famous names =====
\usepackage{xspace}
\newcommand{\Poincare}{Poincar\'e\xspace}
\newcommand{\sch}{Schr\"odinger\xspace}
\newcommand{\SWG}{square waveguide\xspace}%{SWG\xspace}%

%==== colors for editing ======
\newcommand{\change}[1]{{\color{RoyalBlue} #1}}
\newcommand{\comment}[1]{{\color{Maroon} #1}}
\newcommand{\error}[1]{{\color{red} #1}}

% =============================================================================


\begin{document}
\title{Effective optical depth per atom of spin squeezing dynamics using nanophotonic waveguides}
\author{Xiaodong Qi}
\affiliation{Center for Quantum Information and Control, University of New Mexico, Albuquerque, New Mexico 87131, USA}
\author{Ezad Shojaee}
\affiliation{Center for Quantum Information and Control, University of New Mexico, Albuquerque, New Mexico 87131, USA}
\author{Poul S. Jessen}
\affiliation{Center for Quantum Information and Control, University of Arizona, Tucson, Arizona 87521, USA}
\author{Ivan H. Deutsch}
\affiliation{Center for Quantum Information and Control, University of New Mexico, Albuquerque, New Mexico 87131, USA}
\author{Yuan-Yu Jau}
\affiliation{Sandia National Laboratories, Albuquerque, New Mexico 87185, USA}
\date{\today}
\pacs{42.50.Lc, 03.67.Bg, 42.50.Dv, 42.81.Gs}

%================================================================%
\begin{abstract}
We study the strong coupling between photons and atoms that can be achieved in nanophotonic geometries in the dispersive regime to implement efficient quantum interface with neutral atoms.
First, we extend our previous work on spin squeezing and quantum nondemolition (QND) measurement with the birefringence protocol using a nanofiber to a more commonly considered protocol using the Faraday interaction.
We established a general theory to calculate the spin squeezing dynamics and show that $7dB$ of spin squeezing may be achievable with $2500$ atoms trapped $1.8$ times of the fiber radius from the fiber axis, in comparison with $~5dB$ of spin squeezing using the birefringence protocol.
The Faraday interaction protocol does not require a sophisticated search for the magic frequencies and can avoid the difficulty of preparing the atoms on the harder-to-prepare clock states as we did for the birefringence protocol,
and hence is more robust and easier to implement.
Meanwhile, compared to the spin squeezing protocol using Faraday interaction in free space, the nanofiber platform enables us to increase
the photon-atom coupling strength dramatically while to reduce the decoherence due to photon scattering process simultaneously,
and hence can achieve a high peak spin squeezing efficiently.
To achieve an even stronger spin squeezing effect towards a non-Gaussian collective spin state, for example, we propose a generalized optical depth per atom concept applicable to general nanophotonic waveguides and explicitly separate the photon-atom coupling strength from the decoherence characteristic parameter, which can be used to guide the design of novel nanophotonic quantum interfaces and protocols towards efficient quantum control and measurement of atomic states as well as preparing non-Gaussian atomic ensemble states.
We also study the decoherence mechanism for atoms with arbitrary total angular momentum quantum number $f$ with the nanofiber interface, and confirm $f=1$ is the optimal case for the Faraday protocol we study.
Finally, we give an example of analysing a nanophotonic waveguide interface with square cross-section to implement efficiently strong photon-atom coupling and discuss the challenges of using waveguides without cylindrical symmetry beyond optical nanofibers towards efficient photon-atom coupling.
\end{abstract}

\maketitle

%===================INTRODUCTION=====================%
\section{Introduction}

Strong coupling between atoms and photons via nanophotonic structures and such.

The difficulty and challenges.

Our protocol and the structure of this paper.


%========================== Theory ===================================%
\section{Theory} \label{Sec::Theory}

\subsection{Method to calculate spin squeezing dynamics in a QND measurement process}
To study the spin dynamics of the QND measurement process, it is usually valid to calculate the density if the dimension is small. 
For a large ensemble of atoms, however, it requires a big dimension to store the density matrix and it becomes extremely difficult to calculate the evolution of the density operator.
Below, we outline our general method to study the spin dynamics and the spin squeezing evolution without calculating the density operator.

A \textit{squeezed coherent spin state}\index{squeezed spin state!squeezed coherent spin state} is a spin state that the uncertainty principle of two quadratures of collective spin operators saturates and one quadrature is smaller than the other. 
For example, a squeezed coherent spin state may satisfy $ \Delta F_1\Delta F_2=\frac{1}{2} $ (we set $ \hbar=1 $) and $ \Delta F_1<\Delta F_2 $, where $ \hat{F}_1 $ and $ \hat{F}_2 $ are two collective spin operators. 
In general, a \textit{squeezed spin state}\index{squeezed spin state} always satisfies the uncertainty principle relationship $ \Delta F_1\Delta F_2\ge\frac{1}{2} $.
We call the operator that yields the smaller quadrature as the atomic angular momentum operator $ \hat{F}_\perp $ and 
recall the spin squeezing parameter defined by Wineland {\emph{et al.}}~\cite{Wineland1992},
\begin{align}
\zeta^2 &\equiv \frac{\expect{\hat{F}_\parallel(t=0)}^2}{\Delta F_\perp^2(t=0)} \frac{\Delta F_\perp^2}{\expect{\hat{F}_\parallel}^2},
\end{align}
which only requires to calculate the expectation value of the collective atomic angular momentum operator in parallel with the total atomic angular momentum vector in the generalized Bloch sphere, $ \expect{\hat{F}_\parallel} $, as well as the variance of the collective atomic angular momentum operator perpendicular to the total atomic angular momentum operator, $ \Delta F_\perp^2 $. 
Assume the atom number is $ N_A $, these two collective quantities can be decomposed into microscopic quantities by 
\begin{align}
\expect{\Delta F_\perp^2} &= N_A \expect{\Delta f_\perp^2}+\frac{N_A(N_A-1)}{2}\left. \expect{\Delta f_\perp^{(i)}\Delta f_\perp^{(j)}}_s\right|_{i\neq j}\label{eq:DeltaFz2}\\
\expect{\hat{F}_\parallel } &= \sum_i^{N_A} \expect{\hat{f}_\parallel ^{(i)}}=N_A \expect{\hat{f}_\parallel},\label{eq:expectFx}
\end{align}
where the first term of Eq.~\eqref{eq:DeltaFz2} and Eq.~\eqref{eq:expectFx} are solely determined by a symmetric sum over $N_A$ identical spin-$f$ single-body operators, $ \hat{f}_\perp=\hat{f}_\perp^{(i)} $ and $ \hat{f}_\parallel=\hat{f}_\parallel^{(i)} $ with atom labels $ i=1,\cdots,N_A $; the second term of Eq.~\eqref{eq:DeltaFz2} is determined by symmetric two-body covariance terms, $ \left.\expect{\Delta f_\perp^{(i)}\Delta f_\perp^{(j)}}_s\right|_{i\neq j}=\expect{\Delta f_\perp^{(1)}\Delta f_\perp^{(2)}}_s\equiv \expect{\hat{f}_\perp^{(1)}\hat{f}_\perp^{(2)}}_s-\left( \expect{\hat{f}_\perp^{(1)}} \expect{\hat{f}_\perp^{(1)}}\right)_s $, which correspond to the pairwise entanglement among atoms and eventually yield spin squeezing~\cite{Wang2003Spin}.
Above, we have assumed there is a pairwise exchange symmetry among atoms so that we only care about the symmetrized quantities like $ \expect{\Delta f_\perp^{(1)}\Delta f_\perp^{(2)}}_s=\left(\expect{\Delta f_\perp^{(1)}\Delta f_\perp^{(2)}} + \expect{\Delta f_\perp^{(2)}\Delta f_\perp^{(1)}} \right)/2 $. 
Note that the collective state of atoms can be treated as pairwise-symmetric if the detuning is far off resonance or the numbers of atoms is not that large so that the photon scattering among atoms~\cite{Asenjo-Garcia2017Atom,Asenjo-Garcia2017Exponential} can be ignored compared to the measurement backaction which generates spin squeezing as we will discuss later.
In this paper, we will work in the dispersive regime with a few thousands of atoms and ignore the atom-atom interaction caused by photon scattering, and hence the collective atomic system satisfy the exchange symmetry. 

\begin{figure}[htb]
\centering
  \subfloat[b][]{
  \begin{minipage}[bt]{.88\textwidth}
    \centering
    \includegraphics[width=.98\textwidth]{fig/Fig_NanofiberSchematicFaraday}
    \label{fig:Fig_NanofiberSchematicFaraday}
    \end{minipage} }\\
  %\end{subfloat}\\%   
  \subfloat[b][]{
  \begin{minipage}[bt]{.5\textwidth}
      \centering
      \includegraphics[width=.45\textwidth]{fig/poincaresphere_initialS1_crystal}
      \label{fig:poincaresphere_initialS1_crystal}
   \hfill
        \centering
        \includegraphics[width=.45\textwidth]{fig/poincaresphere_initialS1_Faradayrot_crystal}
        \label{fig:poincaresphere_initialS1_Faradayrot_crystal}
   \vfill
    \centering
    \includegraphics[width=.45\textwidth]{fig/blochsphere_initialxFxyz_crystal_orange}
    \label{fig:blochsphere_initialxFxyz_crystal_orange}
   \hfill
      \centering
      \includegraphics[width=.45\textwidth]{fig/blochsphere_initialxFxyz_squeezed_crystal_orange}
      \label{fig:blochsphere_initialxFxyz_squeezed_crystal_orange}
  \end{minipage}}\hfill
  %\end{subfloat}% 
  \subfloat[b][]{
  \begin{minipage}[bt]{.45\textwidth}
        \centering
        \includegraphics[width=.98\textwidth]{fig/stretchedstatestransitions_dashline}
        \label{fig:stretchedstatestransitions}
    \end{minipage}}
    %\end{subfloat}
  \caption{Illustrated are the schematic polarization spectroscopy setup for the QND measurement and spin squeezing generation based on the Faraday effect (a), evolution of the light (top row of (b)) and collective spin states (bottom row of (b)) before and after the measurement, and the internal atomic state transitions in the optical pumping process (c). Diagram (a) is adapted from Ref.~\cite{Qi2016}.}\label{fig:spinsqueezingschematic}
\end{figure}


In the following discussions, we will define the spin squeezing operator $ \hat{F}_\perp=\hat{F}_z=\sum_i\hat{f}_z^{(i)} $. 
We define the fiducial state of an atom as the \textit{up state}, or $ \ket{\uparrow} $. 
Applying $ \hat{f}_z $ on the fiducial state will yield the coupled state or $ \ket{\downarrow} $.
That is, $ \hat{f}_z^{(i)}\ket{\uparrow}^{(i)}\rightarrow \ket{\downarrow}^{(i)} $ for each atom.

To study the spin dynamics, we formally define a stochastic master equation of the atomic ensemble by
\begin{align}\label{eq:totaldrhodt}
\mathrm{d}\hat{\rho}=\left.\mathrm{d}\hat{\rho}\right|_{op} + \left.\mathrm{d}\hat{\rho}\right|_{QND}.
\end{align}
It includes two collective spin dynamic processes. 
The first process is the optical pumping dynamics on each individual atom $i$ positioned at $\br'$ which yields the $\mathrm{d}\hat{\rho}|_{op}=\sum_i^{N_A} \left.\mathrm{d}\hat{\rho}^{(i)}\right|_{op} $ term given by
\begin{align}
&\quad\left.\mathrm{d}\hat{\rho}^{(i)}\right|_{op} =\gamma_s\mathcal{D}^{(i)}\mathrm{d}t\\
&= -\frac{i\gamma_s}{\hbar} \left\{\hat{h}_{\rm eff},\hat{\rho}^{i} \right\}\mathrm{d}t + \gamma_s\sum_q \hat{W}_q(\br')\hat{\rho}^{i}\hat{W}_q^\dagger(\br')\mathrm{d}t,
\end{align}
where the characteristic photon scattering rate $ \gamma_s\equiv \frac{\Gamma_0\Omega^2}{4\Delta_F}=\frac{\sigma_0}{A_{in}}\frac{\Gamma_0^2}{4\Delta_F^2}\dot{N}_L $ with the effective input mode area $ A_{in}=1/n_g|u_{\mathrm{in}}(\br'\!_\perp)|^2 $ and the effective detuning $ \Delta_F $ defined by $ \frac{1}{\Delta_F}=\sum_{f'}\frac{C_{f'ff'}^{(1)}}{\Delta_{ff'}} $ and $ \Delta_{ff'}=\omega-\omega_{ff'} $, where $ \omega_{ff'} $ is the resonance angular frequency between the ground hyperfine structure level $ f $ and the excited hyperfine structure level $ f' $, and $ C_{j'ff'}^{(K)} $ are the coefficients for irreducible rank-$K$ components defined in \cite{Deutsch2010a}.
$\gamma_s$ characterizes the rate of decoherence dynamics and is proportional to the local photon flux of the probe light, $ \dot{N}_L $.

The second term on the right-hand side of Eq.\eqref{eq:totaldrhodt} gives rise to the collective spin dynamics due to QND measurement,
\begin{align}
\left.\mathrm{d}\hat{\rho}\right|_{QND} &= \sqrt{\frac{\kappa}{4}}\mathcal{H}\left[\hat{\rho} \right]\mathrm{d}W + \frac{\kappa}{4}\mathcal{L}\left[ \hat{\rho}\right]\mathrm{d}t.
\end{align}
Above, we have defined the measurement strength $\kappa \equiv |\chi|^2\dot{N}_L\equiv \frac{\sigma_0A_{in}}{A_{int}^2}\gamma_s $ determining the rate of the spin squeezing in absence of decoherent processes, where $\dot{N}_L$ is the photon number flux, $\chi$ is the light-atom coupling strength and $A_{int}$ is the effective interaction mode area which can be specified for a particular QND measurement protocol. We have also assumed the measurement backation is a stochastic Weiner process where $\mathrm{d}W$ is the increment satisfying $\mathrm{d}W^2 = \mathrm{d}t$. The conditional dynamics responding to the measurement evolve under the superoperator
\begin{align}
\mathcal{H}\left[ \hat{\rho}\right] &= \hat{F}_\perp\hat{\rho} + \hat{\rho}\hat{F}_\perp -2\expect{\hat{F}_\perp}\hat{\rho}
\end{align}
and the collective Lindblad map due to the direct photon scattering of the guided modes from the atoms
\begin{align}
\mathcal{L}\left[ \hat{\rho} \right] &= \hat{F}_\perp\hat{\rho}\hat{F}_\perp-\frac{1}{2}\left(\hat{\rho}\hat{F}_\perp^2+\hat{F}_\perp^2\hat{\rho} \right)=\frac{1}{2}\left[\hat{F}_\perp,\left[\hat{\rho},\hat{F}_\perp \right] \right].
\end{align}

As shown in the equations above, the spin squeezing dynamics is a competition between the coherent squeezing process and all decoherent processes which are characterized by $\kappa$ and $\gamma_s$, respectively. 
If we define an effective cooperativity or optical depth (OD) per atom quantity for the spin squeezing dynamics by
\begin{align}
\frac{\mathrm{OD}_{\rm eff}}{N_A} \equiv \frac{\kappa}{\gamma_s}=\frac{\sigma_0A_{in}}{A_F^2},
\end{align}
the peaking spin squeezing dynamics can then be characterized by $\frac{\mathrm{OD}_{\rm eff}}{N_A}$, and the geometry of the spin squeezing protocol can then be roughly designed with the goal to maximize $\frac{\mathrm{OD}_{\rm eff}}{N_A}$ by minimizing $A_{in}$ and maximizing $A_{int}$.  

The full expression for the cooperativity is
\begin{align}\label{eq:c1_bound}
C_1(\mbf{r}'_\perp) = n_g\sigma_0\frac{ \left\vert \text{Re} \left[ \mbf{u}^*_V (\mbf{r}'_\perp) \times \mbf{u}_H (\mbf{r}'_\perp) \right]\right\vert^2}{\vert \mbf{u}_H (\mbf{r}'_\perp) \vert ^2} \le n_g  \sigma_0 \vert \mbf{u}_V (\mbf{r}'_\perp) \vert ^2.
\end{align}
The upper bound is typically not achieved for a nanophotonic waveguide due to the out-of-phase $z $-component of the modes at arbitrary positions.  In free-space, or for approximate TEM modes of a waveguide, $ C_1\propto\sigma_0/\Ai=\sigma_0 u(r'_\perp)^2$, where $ \vert \mbf{u}_H (\mbf{r}'_\perp) \vert ^2=\vert \mbf{u}_V (\mbf{r}'_\perp) \vert ^2 = u(r'_\perp)^2 $ independent of azimuthal angle. 
%The upper bound may also be saturated with pure TEM modes of a waveguide when the longitudinal mode component vanishes. 
Even with the longitudinal component, however, by proper choice of the atomic position for the waveguide geometry, the cooperativity can be strongly enhanced over free space, due to the tight and anisotropic confinement of the modes the waveguide guides,  as we  show below.

We can bring in the spin coherence operator $\hat{\sigma}_{ba}=\ket{b}\bra{a}$ to represent the matrix element of any atomic angular momentum operator $ \hat{f}_m $ ($ m=x,y,z $) of a single atom, and hence the spin squeezing dynamics can be characterized by the expectation value of single-body spin coherence operators $\expect{\hat{\sigma}_{ba}}$ and the symmetric two-body covariances $\expect{\Delta \sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)} }_s$. 
%the symmetric three-body correlations $\expect{\Delta \sigma^{(1)}_{b_1a_1}\Delta \sigma^{(2)}_{b_2a_2}\Delta \sigma^{(3)}_{b_3a_3} }_s$ and so on. 
%As higher-order correlations becomes negligible, one can cut off the correlation terms at a certain order.
If one can truncate the spin dynamics up to the two-body correlations, we only need the following two sets of stochastic differential equations:
\begin{subequations}
\begin{align}
d\expect{\hat{\sigma}_{ba}} &=\left. d{\expect{\hat{\sigma}_{ba}}}\right|_{op} + \left. d{\expect{\hat{\sigma}_{ba}}}\right|_{\mathcal{H}}+\left. d{\expect{\hat{\sigma}_{ba}}}\right|_{\mathcal{L}} \\
d\expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{dc}^{(2)}}_s &= \left. d{\expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{dc}^{(2)}}_s}\right|_{op} + \left. d{\expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{dc}^{(2)}}_s}\right|_{\mathcal{H}} + \left. d{\expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{dc}^{(2)}}_s}\right|_{\mathcal{L}}.
\end{align}
\end{subequations}

In details, the optical dynamics part can be given by
\begin{align}
\left. \dt{\expect{\hat{\sigma}_{ba}}}\right|_{op} &= \gamma_s\expect{\mathcal{D}^\dagger \left[ \hat{\sigma}_{ba}\right]}\\
&= \gamma_s\sum_{d,c}\tr\left(\mathcal{D}^\dagger \left[ \hat{\sigma}_{ba}\right]\hat{\sigma}_{dc} \right)\expect{\hat{\sigma}_{dc} }\\
\left. \dt{\expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{dc}^{(2)}}_s}\right|_{op} &=\gamma_s\expect{\Delta\mathcal{D}^\dagger[\hat{\sigma}_{ba}^{(1)}]\Delta\sigma_{dc}^{(2)} }_s + \gamma_s\expect{\Delta\sigma_{ba}^{(1)}\Delta\mathcal{D}^\dagger[\hat{\sigma}_{dc}^{(2)}] }_s\\
&= \gamma_s\sum_{m,n}\tr\left(\mathcal{D}^\dagger[\hat{\sigma}_{ba}]\hat{\sigma}_{mn} \right)\expect{\Delta \sigma_{mn}^{(1)}\Delta \sigma_{dc}^{(2)} }_s + \gamma_s\sum_{m,n}\tr\left(\mathcal{D}^\dagger[\hat{\sigma}_{dc}]\hat{\sigma}_{mn} \right) \expect{\Delta \sigma_{ba}^{(1)}\Delta \sigma_{mn}^{(2)} }_s.
\end{align} 
Similarly, we will need the one- and two-body correlations due to the $ \mathcal{H} $ and $ \mathcal{L} $ superoperators given by the following.
\begin{subequations}
\begin{align}
\left.d\expect{\hat{\sigma}_{ba}}\right|_\mathcal{H} &=\sqrt{\frac{\kappa}{4}}\expect{\mathcal{H}^\dagger\left[\hat{\sigma}_{ba} \right]}dW \\
\left.d\expect{\hat{\sigma}_{ba}}\right|_\mathcal{L} &= \frac{\kappa}{4}\expect{\mathcal{L}^\dagger\left[\hat{\sigma}_{ba} \right]}dt
\end{align}
\end{subequations}
In principle, the two-body covariance terms can be coupled to high-order many-body terms. 
In our case, we assume the state of the ensemble can be well captured in the symmetric Gaussian state limit, and hence the two-body covariance equations due to the collective measurement can be given by
\begin{subequations}
\begin{align}
\left.d\expect{\Delta \sigma_{ba}^{(1)} \Delta \sigma_{dc}^{(2)}}_s \right|_\mathcal{H} &= -\kappa\expect{\Delta\sigma_{ba}^{(1)}\Delta F_\perp }_s \expect{\Delta F_\perp \Delta \sigma_{dc}^{(2)} }_s dt \\
\left.d\expect{\Delta \sigma_{ba}^{(1)} \Delta \sigma_{dc}^{(2)}}_s\right|_\mathcal{L} &= 0.
\end{align}
\end{subequations}

\subsection{A spin squeezing protocol using Faraday interaction with spin coherent states}

\begin{figure}[htb]
\centering
  \includegraphics[width=.69\textwidth]{fig/nanofiberswg_Hmode6} 
  \caption{$H$-mode components and intensity distribution in the $ xy $-plane for the optical nanofiber and SWG. Shown on the first row, from left to right, are the real part of $ E_x(\br\!_\perp) $, the real part of $ E_z(\br\!_\perp) $ and the imaginary part of $ E_z(\br\!_\perp) $ for the nanofiber. White lines outline the waveguide boundary. The scale is normalized to the maximum value of all field components shown on the subplots. Plots on the second row are the same components as the first row but for the SWG. All the other parts of the electric field components of the $ H $- or TE$_{01}$-mode not shown in these plots for both waveguide geometries are zero everywhere. On the third row, we plot the normalized intensity distribution on the transverse plane for both geometries. The arrows indicates the local field's direction and amplitude (relative length) along the vertical waveguide axis, which only have an $ x $-component. The $ V $-mode is the $ H $-mode rotated by $ 90^\circ $ around the waveguide axis. }\label{fig:nanofiberswg_E_ints}
\end{figure}


Following the process demonstrated in our previous work~\cite{Qi2016}, the light-atom interaction Hamiltonian with one atom can be written as
\begin{align}
\hat{h}_\eff &= -\hat{\mathbf{E}}^{(-)}(\br')\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot\hat{\mathbf{E}}^{(+)}(\br')\nn\\
&= -\frac{2\pi\hbar\omega}{v_g}\left[\mathbf{u}_H^*\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot \mathbf{u}_H\hat{a}_H^\dagger\hat{a}_H\right.
+ \mathbf{u}_H^*\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot \mathbf{u}_V\hat{a}_H^\dagger\hat{a}_V\nn\\
&\quad\quad + \mathbf{u}_V^*\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot \mathbf{u}_H\hat{a}_V^\dagger\hat{a}_H 
\left. + \mathbf{u}_V^*\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot \mathbf{u}_V\hat{a}_V^\dagger\hat{a}_V\right]\\
&= \hbar\left[(\hat{\chi}_{HH}+\hat{\chi}_{VV})\hat{S}_0 + (\hat{\chi}_{HH}-\hat{\chi}_{VV})\hat{S}_1 + (\hat{\chi}_{HV}+\hat{\chi}_{VH})\hat{S}_2 + i(\hat{\chi}_{HV}-\hat{\chi}_{VH})\hat{S}_3 \right]\\
%\hbar \left[\left(\chi_{RR\uparrow} + \chi_{RR\downarrow} +\chi_{LL\uparrow}+\chi_{LL\downarrow} \right)\hat{F}_0\hat{S}_0 \right.\nonumber\\
%&\quad+\left(\chi_{RR\uparrow} + \chi_{RR\downarrow} -\chi_{LL\uparrow}-\chi_{LL\downarrow} \right)\hat{F}_0\hat{S}_3\nonumber\\
%&\quad+\left(\chi_{RR\uparrow} + \chi_{LL\uparrow} -\chi_{RR\downarrow}-\chi_{LL\downarrow} \right)\hat{F}_3\hat{S}_0\nonumber\\
%&\quad+\left(\chi_{RR\uparrow} - \chi_{RR\downarrow} +\chi_{LL\downarrow}-\chi_{LL\uparrow} \right)\hat{F}_3\hat{S}_3\nonumber\\
%&\quad+i\left(\chi_{LR\uparrow} - \chi_{RL\uparrow} +\chi_{RL\downarrow}-\chi_{RL\downarrow} \right)\hat{F}_0\hat{S}_1\nonumber\\
%&\quad+\left(\chi_{RL\uparrow} + \chi_{LR\uparrow} +\chi_{RL\downarrow}+\chi_{LR\downarrow} \right)\hat{F}_0\hat{S}_2\nonumber\\
%&\quad+i\left(\chi_{LR\uparrow} - \chi_{RL\uparrow} +\chi_{RL\downarrow}-\chi_{LR\downarrow} \right)\hat{F}_3\hat{S}_1\nonumber\\
%&\quad+\left.\left(\chi_{LR\uparrow} + \chi_{RL\uparrow} -\chi_{LR\downarrow}-\chi_{RL\downarrow} \right)\hat{F}_3\hat{S}_2 \right]\\
&=\hbar\sum_{i=0}^3 \hat{\chi}_{i}\hat{S}_i\\
&=\hbar\sum_{i,j=0} \chi_{ij}\hat{f}_i\hat{S}_j,
\end{align}
where $ \hat{S}_i $ are the Stokes vector operators of the light indicating its polarization, and the mode-atom coupling operator
\begin{align}
\hat{\chi}_{pp'} 
&=-\frac{2\pi \omega}{v_g}\mathbf{u}_{p}^*(r'\!_\perp,\phi')\cdot \hat{\tensor{\alpha}}\cdot \mathbf{u}_{p'}(r'\!_\perp,\phi')\\
&= \sum_{f'} \frac{n_g\sigma_0}{4}\frac{\Gamma_{f'}}{\Delta_{ff'}+i\Gamma_{f'}/2}\cdot \left\{ C_{j'ff'}^{(0)}\mathbf{u}_p^*(r'\!_\perp)\cdot \mathbf{u}_{p'}(r'\!_\perp)\hat{\mathbbm{1}}\right.\nn\\
&\quad\quad +iC_{j'ff'}^{(1)}\left(\mathbf{u}_p^*(r'\!_\perp)\times\mathbf{u}_{p'}(r'\!_\perp) \right)\cdot \hat{\mathbf{f}} \nonumber\\
&\quad\quad\left. + C_{j'ff'}^{(2)}\sum_{i,j}\left[u^*_{p,i}u_{p',j}(\frac{\hat{f}_i\hat{f}_j+\hat{f}_j\hat{f}_i}{2}-\frac{\delta_{ij}}{3}\hat{\mathbf{f}}\cdot\hat{\mathbf{f}}) \right]\right\}
%&\left.+C_{jj'ff'}^{(2)}\left[\mathbf{u}_p^*(r'\!_\perp)\cdot \mathbf{u}_{p'}(r'\!_\perp)\left(\frac{f(f+1)}{6}-\frac{m^2}{2} \right)+\mathbf{u}_p^*(r'\!_\perp)\cdot (\hat{e}^*_{\tilde{z}}\hat{e}_{\tilde{z}})\cdot \mathbf{u}_{p'}(r'\!_\perp)\left(\frac{3m^2}{2}-\frac{f(f+1)}{2} \right) \right] \right\}
\label{eq:chippp}
\end{align}
with the horizontally(H)- and vertically(V)-linearly polarized guided modes, $ \mathbf{u}_p(r'\!_\perp) $, at the atom position $ \br'=(r'\!_\perp,\phi',z') $. 
$ \chi_{ij}=\tr[\hat{f}_i\hat{\chi}_j]/(2f+1) $ is the coupling strength between spin operator $ \hat{f}_i $ and Stokes operator $ \hat{S}_j $. 
For example, $ \chi_{33} $ is the coupling strength between $ \hat{f}_z $ and $ \hat{S}_3 $.
The fundamental guided modes of an optical nanofiber has been defined in the appendix of our previous paper~\cite{Qi2016}. 
In general, for a cylindrical waveguide, the H- and V-modes are the guided modes adiabatically transferred from a corresponding linearly polarized input light from one end of the waveguide, where H- and V-directions are orthogonal to each other in the transverse plane.
The coupling operator or Eq.\eqref{eq:chippp} includes three terms corresponding to scalar, vector and tensor interactions between atoms and the probe light which are proportional to $ C_{j'ff'}^{(K)} $ with $ K=0,\,1,\,2 $, respectively.

Faraday interaction is the interaction when the helicity of the input and output of the light signal is preserved. 
In the context of QND measurement, a Faraday interaction based protocol is the protocol when the phase change of the output light on the equator of the \Poincare sphere with a linear polarization input is used to calibrate the spin state after the atom-light interaction.
This phase change of light polarization corresponds to a rotation about $ \mathbf{S}_3 $ axis on the \Poincare sphere. 
From Eq.\eqref{eq:chippp}, the coupling atomic term is proportional to 
\begin{align}
\hat{\chi}_{i3} &= \hat{\chi}_{HV}-\hat{\chi}_{VH}=\hat{\chi}_{RR}-\hat{\chi}_{LL},
\end{align}
where the last step is derived in appendix~\ref{Appendix:LRbases}. These relationships indicate the Faraday effect is generated by the intensity or photon number difference between the right- and left-polarized modes or by the phase difference between two linearly polarized modes due to the polarizability of the atoms.
One can prove that the strength of the Faraday interaction is dominated by the vector interaction term in Eq.\eqref{eq:chippp}.
This implies that, to maximize the Faraday interaction, it is optimal to choose a quantization axis along the direction of $ \mathbf{v}_F=\mathbf{u}_H^*(\br'\!_\perp)\!\times\!\mathbf{u}_{V}(\br'\!_\perp)\!-\!\mathbf{u}_V^*(\br'\!_\perp)\!\times\!\mathbf{u}_{H}(\br'\!_\perp)$ given the atoms are placed at $ \br'\!_\perp $ position in the transverse plane of the waveguide. 
In reality, this product of modes may be elliptical with at least one direction component is imaginary while others are real.
If this happens, since the quantization axis is the direction in which a magnetic field is pointing to in 3D real space, the optimal choice of quantization axis should be the direction corresponding to the largest component of the $ \mathbf{v}_F $ vector. 
For a cylindrical waveguide, this doesn't seem to happen.
Take the example of a nanofiber, at an arbitrary position $ \br'\!_\perp=(r'\!_\perp,\phi') $ of atoms in the transverse plane,
\begin{align}
\mathbf{u}_H^*(r'\!_\perp,\phi')\times \mathbf{u}_V(r'\!_\perp,\phi') &= 2u_{r\!_\perp} u_\phi\mathbf{e}_z - 2iu_zu_{r\!_\perp}\sin2\phi \mathbf{e}_\phi + 2iu_\phi u_z\cos2\phi \mathbf{e}_{r\!_\perp} \\
\mathbf{u}_V^*(r'\!_\perp,\phi')\times \mathbf{u}_H(r'\!_\perp,\phi') &= -2u_{r\!_\perp} u_\phi\mathbf{e}_z - 2iu_zu_{r\!_\perp}\sin2\phi \mathbf{e}_\phi + 2iu_\phi u_z\cos2\phi \mathbf{e}_{r\!_\perp},
\end{align}
and therefore,
\begin{align}\label{eq:Faradayaxis}
\mathbf{v}_F=\mathbf{u}_H^*(\br'\!_\perp)\!\times\!\mathbf{u}_{V}(\br'\!_\perp)\!-\!\mathbf{u}_V^*(\br'\!_\perp)\!\times\!\mathbf{u}_{H}(\br'\!_\perp) = 4u_{r\!_\perp} u_\phi\mathbf{e}_z,
\end{align}
where $ u_{r\!_\perp}=u_{r\!_\perp}(r'\!_\perp) $, $ u_\phi=u_\phi(r'\!_\perp) $ and $ u_z=u_z(r'\!_\perp) $ are the right-circularly polarized mode components independent of longitudinal and azimuthal positions defined in the appendix A of Ref~\cite{Qi2016} or the corresponding appendix in this dissertation.
Based on this result, choosing $ z $-direction as the quantization axis is optimal for the QND measurement and spin squeezing protocol for atoms trapped near an optical nanofiber.

This conclusion can be generalized to cylindrical waveguides at large, which have a smooth or slow change of index of refraction along the light propagation direction in the wavelength scale while the cross-section of the waveguides can be arbitrary.
From the perspective of transformation optics~\cite{Leonhardt2006Optical,Kundtz2011Electromagnetic}, we can consider a set of $H$- and $V$-modes for the new waveguide are generated by adiabatically transforming the orthogonal set of modes from a cylindrical nanofiber to the target cross-section shape of the waveguide.
Since the transformation is approximately limited to the $xy$-plane of the coordinate system transformation determined by a Jacobian matrix, Eq.\eqref{eq:Faradayaxis} will preserve the form in the new waveguide coordinate system where only $ z $-component is non-zero and should be chosen as the optimal choice of the quantization axis.

With the quantization axis chosen along $ \mathbf{e}_{\tilde{z}}=\mathbf{e}_z $, one can define the effective QND measurement Hamiltonian for an atomic ensemble by
\begin{align}
\hat{H}_F &= \hbar \chi_{33}\hat{F}_z \hat{S}_3,
\end{align}
where $ \chi_{33} $ is the measurement strength characterizing the entanglement between the collective spin state of $ \hat{F}_z $ and the probe's polarization state of $ \hat{S}_3 $.
Now that $ \hat{F}_\perp=\hat{F}_z $ is the squeezing operator. 
Let us consider an ensemble of $ ^{133} $Cs atoms initially prepared as a spin coherent state (SCS) with every atom in the stretch state of the $ 6S_{1/2}$ $f=4 $ ground manifold in the $ x $-basis, where the quantization $ x $-direction is along the diagonal or $ \phi=\pi/4 $ direction in the $ H $-$ V $ Cartesian coordinate system sitting on the fiber axis.
One can show that a stretch state can generate the maximum coupling between atoms and light when all atoms are prepared in a SCS.
We call this initial state of one atom as the fiducial state or $ \ket{\phi_0}=\ket{\uparrow}_x = \ket{6S_{1/2},f=4,m_x=4} $ and the collective SCS can be written as $ \ket{\Psi_0}=\ket{\uparrow}_x^{\otimes N_A} $.
We define the coupled state by applying the individual squeezing operator $ \hat{f}_z$ on the fiducial state $\ket{\uparrow}_x $.
In our case, the coupled state can be defined as $ \ket{\downarrow}_x=\ket{6S_{1/2},f=4,m_x=3} $.
Similarly, to include the transfer of coherence in the squeezing process, we define the transfer state as $ \ket{T}_x=\ket{6S_{1/2},f=4,m_x=2} $ resulted from applying $ \hat{f}_z$ on the coupled state. 
For simplicity, we remove all the subscript $ x $ of quantum states and assume we always work in the $ x $-basis if no explicit notations for the Faraday interaction protocol. 

We consider the far-detuning regime, where we may be able to set the decay rates of the excited levels $ \Gamma_{f'}= \Gamma_0$ as constant for all $ f' $ in the same fine structure manifold, where we consider $ \Gamma_0 $ as the averaged modified decay rates from the excite fine structure level to the ground level when atoms are in a completely mixed state for simplicity.
We also ignore the tensor coupling strength related to $ C_{jj'ff'}^{(2)} $ terms in Eq.\eqref{eq:chippp} as the tensor interaction strength ($ \sim 1/\Delta^2 $) is relatively small compared to the vector interaction strength ($ \sim 1/\Delta $)~\cite{Deutsch2010a}. 
For a nanofiber geometry, the Faraday interaction coupling strength is independent of the azimuthal position of the atoms and can be simplified as
\begin{align}
\chi_{33} &= -\sum_{f'}n_g\sigma_0\frac{\Gamma_0}{\Delta_{ff'}+i\Gamma_0/2}C_{jj'ff'}^{(1)}u_{r\!_\perp}(r'\!_\perp)u_\phi(r'\!_\perp)\\
&=\frac{\sigma_0}{A_F}\frac{\Gamma_0}{\Delta_F},
\end{align}
where the effective Faraday interaction mode area $ A_F=1/2n_g|u_{r\!_\perp}(r'\!_\perp)u_\phi(r\!_\perp)| $, and the effective detuning $ \Delta_F=\sum_{f'}\frac{-C_{j'ff'}^{(1)}}{\Delta_{ff'}} $.
The measurement strength is now defined as
\begin{align}
\kappa\equiv|\chi_{33}|^2\dot{N}_L=\frac{\sigma_0A_{in}}{A_F^2}\gamma_s,
\end{align}
where the characteristic photon scattering rate $ \gamma_s\equiv \frac{\Gamma_0\Omega^2}{4\Delta_F}=\frac{\sigma_0}{A_{in}}\frac{\Gamma_0^2}{4\Delta_F^2}\dot{N}_L $ and the effective mode area $ A_{in}=1/n_g|u_{\mathrm{in}}(\br'\!_\perp)|^2 $.
Now we can define the OD per atom for the Faraday interaction using SCS by
\begin{align}
\frac{\mathrm{OD}}{N_A} \equiv \frac{\kappa}{\gamma_s}=\frac{\sigma_0A_{in}}{A_F^2}.
\end{align}

Since the Faraday measurement strength doesn't depend on the azimuthal direction of the atom position, ideally, to implement an optimal Faraday interaction geometry, it is preferable to place atoms along some azimuthal direction with the minimum impact from birefringence effects as well as the decoherence due to the presence of the nanofiber. 

Firstly, in the far-detuning regime, since the birefringence effect is dominated by the scalar coupling which is proportional to the intensity difference between the $H$ ($x$)- and $V$ ($ y $)-mode components at the atom positions, both diagonal and anti-diagonal directions in the transverse plane yields a minimum birefringence effect where the intensity of the $ H $- and $ V $-mode components are equal, based on the symmetry of the fiber and given a diagonally polarized $ D $-mode input.
That is the optimal position of atoms with vanished birefringence effect due to the intensity difference of the local $ H $- and $ V $-modes could be along the $ \phi'=n\pi/4 $ ($ n=1,3,5,7 $) radial directions. 
Secondly, to minimize the decoherence damages to spin squeezing, the optimal position of the atom should be chosen so that $ A_{\mathrm{in}} $ reaches the minimum value, which yields $ \phi'=3\pi/4 $ or $ 7\pi/4 $ where $ A_{\mathrm{in}}=1/2n_g|u_\phi(r'\!_\perp)|^2 $.
Combining these two factors, we find the optimal choice of atoms' azimuthal positions are along the anti-diagonal direction, that is $\phi'=3\pi/4 $ or $ 7\pi/4 $.

The number of equations needed to include different subspace...

\subsection{Cancel tensor light-shift with two-color probes}



%=================== Simulations with waveguides =====================%
\section{Discussions} \label{Sec::Discussions}
\comment{Leave discussions on how does the internal structure of atoms and the total atomic quantum number affect the spin squeezing parameter for future works.}

% set package options
\captionsetup[subfloat]{position=top,singlelinecheck=off,labelfont={normalsize,sf}, %justification=raggedright,
  labelformat=simple,listofformat=subparens,aboveskip=0pt,parskip=0pt,farskip=5pt,
  captionskip=0pt}

% customize subfigure label to capitals
\renewcommand{\thesubfigure}{(\alph{subfigure})}

%\newlength\figureheight
%\newlength\figurewidth
\begin{figure}[htb]
\centering
%\pgfplotsset{yticklabel style={text width=3em,align=right}}
 \begin{minipage}[h]{0.2\linewidth}
 %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
  \subfloat[h][$1/\AF$]{
    %\setlength\figureheight{0.3\textwidth}
    %\setlength\figurewidth{0.3\textwidth}
    \input{fig/nanofiber_invA_Far_xy.tex}\label{fig:nanofiber_invAFarxy}
    }
    \hfill
  \subfloat[h][$ 1/\Ain $]{
      %\setlength\figureheight{0.3\textwidth}
      %\setlength\figurewidth{0.3\textwidth}
      \label{fig:nanofiber_invAin_xy}
      \input{fig/nanofiber_invA_in_xy.tex}
      }
   \end{minipage}\hfill
   \begin{minipage}[h]{0.75\linewidth}
   \subfloat[h][$ C_1 $]{
      %\setlength\figureheight{0.6\textwidth}
      %\setlength\figurewidth{0.6\textwidth}
      \label{fig:nanofiber_C1_xy}
      \input{fig/nanofiber_C1_xy.tex}
      }
   \end{minipage}
   %\end{tabular}
\caption{Distributions of effective mode areas and cooperativity per atom near an optical nanofiber.}\label{fig:nanofiber_Aeffgeometry}
\end{figure}

\begin{figure}[htb]
\centering
%\pgfplotsset{yticklabel style={text width=3em,align=right}}
 \begin{minipage}[h]{0.2\linewidth}
 %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
  \subfloat[h][$1/\AF$]{
    %\setlength\figureheight{0.3\textwidth}
    %\setlength\figurewidth{0.3\textwidth}
    \input{fig/swg_invA_Far_xy.tex}\label{fig:swg_invAFarxy}
    }
    \hfill
  \subfloat[h][$ 1/\Ain $]{
      %\setlength\figureheight{0.3\textwidth}
      %\setlength\figurewidth{0.3\textwidth}
      \label{fig:swg_invAin_xy}
      \input{fig/swg_invA_in_xy.tex}
      }
   \end{minipage}\hfill
   \begin{minipage}[h]{0.75\linewidth}
   \subfloat[h][$ C_1 $]{
      %\setlength\figureheight{0.6\textwidth}
      %\setlength\figurewidth{0.6\textwidth}
      \label{fig:swg_C1_xy}
      \input{fig/swg_C1_xy.tex}
      }
   \end{minipage}
   %\end{tabular}
\caption{Distributions of effective mode areas and cooperativity per atom near a SWG.}\label{fig:swg_Aeffgeometry}
\end{figure}

%=================== A toy model with a spin-1/2 system =====================%
\subsection{Spin squeezing with a spin-$1/2$ system} \label{Sec::squeezingwithspinhalfsystems}


%============================ spin squeezing with nanofibers ==============================%
\subsection{QND measurement and spin squeezing with nanofibers} \label{Sec::Nanofiber}

The optimal choice of atom position.

Given an ensemble of $N_A$ atoms initially prepared in a spin coherent state for the hyperfine spin $f$, polarized in the transverse plane, e.g., along the $x$-axis, a QND measurement of the collective spin  $F_z$ will squeeze the uncertainty of that component.  The metrologically relevant squeezing parameter~\cite{Wineland1992} is,
\begin{align}\label{eq:xi2Faraday}
\xi^2 &\equiv  \frac{2 N_A f \expect{\Delta F_z ^2}}{\expect{\hat{F}_x}^2}.
\end{align}
Under the assumption that the state is symmetric with respect to exchange of any two atoms, valid when we start in a spin coherent state all coupling is uniform over the ensemble, the collective expectation value can be decomposed into 
\begin{subequations}\label{eq:Ftof_squeezing}
\begin{align}
\expect{\Delta F_z^2} &= N_A \expect{\Delta f_z^2}+N_A(N_A-1)\left. \expect{\Delta f_z^{(i)}\Delta f_z^{(j)}}\right|_{i\neq j}\label{eq:DeltaFz2}\\
\expect{\hat{F}_x } & =N_A \expect{\hat{f}_x} \label{eq:expectFx}.
\end{align}
\end{subequations}
 The first term of Eq.~\eqref{eq:DeltaFz2} and Eq.~\eqref{eq:expectFx} is the projection noise associated with the  $N_A$ identical spin-$f$  atoms, and  the second term of Eq.~\eqref{eq:DeltaFz2} is determined by two-body covariances, $ \left.\expect{\Delta f_z^{(i)}\Delta f_z^{(j)}}\right|_{i\neq j}=\expect{\Delta f_z^{(1)}\Delta f_z^{(2)}} = \expect{\hat{f}_z^{(1)}\hat{f}_z^{(2)}}-\expect{\hat{f}_z^{(1)}} \expect{\hat{f}_z^{(1)}} $.  Negative values in these two-body correlations correspond to the pairwise entanglement among atoms that yields spin squeezing~\cite{Wang2003Spin}.  Note, when detuning is sufficiently far-off resonance so that all collective sub- and super-radiant modes~\cite{Asenjo-Garcia2017Atom,Asenjo-Garcia2017Exponential} are equally (and thus symmetrically) excited.  In this paper, we work in the dispersive regime with a few thousands of atoms and can safely ignore the atom-atom interaction caused by multiple scattering, and hence the collective atomic system satisfy the exchange symmetry. 

To study the spin squeezing dynamics, we employ a first-principles stochastic master equation for the collective state of $N_A$ atoms,
\begin{align}\label{eq:totaldrhodt}
\mathrm{d}\hat{\rho}= \left.\mathrm{d}\hat{\rho}\right|_{QND}+\left.\mathrm{d}\hat{\rho}\right|_{op}.
\end{align}
The first term on the right-hand side of Eq.\eqref{eq:totaldrhodt} governs the spin dynamics arising from QND measurement~\cite{Jacobs2006,Baragiola2014},
\begin{align}
\left.\mathrm{d}\hat{\rho}\right|_{QND} &= \sqrt{\frac{\kappa}{4}}\mathcal{H}\left[\hat{\rho} \right]\mathrm{d}W + \frac{\kappa}{4}\mathcal{L}\left[ \hat{\rho}\right]\mathrm{d}t, 
\end{align}
where  $\kappa$ is the measurement strength defined in Eq.~\eqref{eq:kappa}, and $\mathrm{d}W$ is a stochastic Weiner interval. The conditional dynamics are generated by superoperators that depend on the {\em collective} spin
\begin{subequations}
\begin{align}
\mathcal{H}\left[ \hat{\rho}\right] &= \hat{F}_z \hat{\rho} + \hat{\rho}\hat{F}_z -2\expect{\hat{F}_z}\hat{\rho}, \\
\mathcal{L}\left[ \hat{\rho} \right] &= \hat{F}_z \hat{\rho}\hat{F}_z -\frac{1}{2}\left(\hat{\rho}\hat{F}_z^2+\hat{F}_z^2\hat{\rho} \right)=\frac{1}{2}\left[\hat{F}_z,\left[\hat{\rho},\hat{F}_z \right] \right].
\end{align}
\end{subequations}
The second term governs decoherence arising from optical pumping, which acts {\em locally} on each atom$,\mathrm{d}\hat{\rho}|_{op}=\sum_n^{N_A} \mathcal{D}^{(n)}\left[ \hat{\rho}\right] \mathrm{d}t$, where 
\begin{equation}
\mathcal{D}^{(n)}\left[ \hat{\rho}\right] = -\frac{i}{\hbar}\left(\hat{H}^{(n)}_{\rm eff}\hat{\rho} - \hat{\rho} \hat{H}^{(n)\dag}_{\rm eff}\right) + \gamma_{op} \sum_q \hat{W}^{(n)}_q \hat{\rho}\hat{W}^{(n)\dag}_q.
\label{op_superator}
\end{equation}
Here $\hat{H}^{(n)}_{\rm eff}$ is the effective nonHermitian Hamiltonian describing the local light shift and absorption by the $i^{th}$ atom and $\hat{W}^{(n)}_q$ is the jump operator corresponding to optical pumping through spontaneous emission of a photon of polarization $q$~\cite{Deutsch2010a} (see Appendix~\ref{Sec::opticalpumpinginrotatingframe}) \comment{(we never defined the jump operators in this paper. These jump operators might have a different form from Ref.\cite{Deutsch2010a} as we are using $ \gamma_{op} $ here.)}.   
The rate of decoherence is characterized by optical pumping rate, $\gamma_{op}$.  Note, the optical pumping superoperator, Eq.\eqref{op_superator},  is not trace preserving when restricted to a given ground-state hyperfine manifold $f$.  In this case, optical pumping of atoms to the other hyperfine manifold in the ground-electronic state is treated as loss.  If the atoms are placed at the optimal position, the local field is linearly polarized.  In that case the vector light shift vanishes, and for detunings large compared to the excited-state hyperfine splitting, the rank-2 tensor light shift is negligible over the time scales of interest.  In that case the light shift is dominated by the scalar component, which has no effect on the spin dynamics.  In that case $\hat{H}_{\rm eff} = -i\hbar \gamma_{op}/2 1$.

The solution to the master equation is made possible by three approximations. Firstly, we restrict the subspace of internal magnetic sublevels that participate in the dynamics.  The system is initialized in a spin coherent state, with all atoms spin-polarized along the $x$-axis.  We denote this as the ``fiducial state" $\ket{\uparrow} = \ket{f, m_x =f}$.   Through QND measurement, spin squeezing is induced by entanglement with the  ``coupled state"  $\ket{\downarrow} = \ket{f, m_x=f-1}$.  Optical pumping is dominated by ``spin flips" $\ket{\uparrow}\rightarrow \ket{\downarrow}$ and ``loss" due to pumping to the other hyperfine level.  We finally include a third internal magnetic sublevel $\ket{T} = \ket{f, m_x=f-2}$ to account for  ``transfer of coherences" that can occur in spontaneous emission~\cite{Norris2012Enhanced,Norris2014} (see Fig.~\ref{fig:spinsqueezinglevelstructure}).  Restricted to this qutrit basis, the internal hyperfine spin operators are
\begin{subequations}\label{eq:f_in_xbasis}
\begin{align}
\hat{f}_x &= f \hat{\sigma}_{\uparrow \uparrow} +(f-1) \hat{\sigma}_{\downarrow \downarrow} + (f-2)  \hat{\sigma}_{T T}, \\
\hat{f}_y &=-i \sqrt{\frac{f}{2}} \left(\hat{\sigma}_{\uparrow \downarrow} - \hat{\sigma}_{\downarrow \uparrow}\right) -i \sqrt{\frac{2f-1}{2}}  \left(\hat{\sigma}_{\downarrow T} - \hat{\sigma}_{T \downarrow }\right) \\
\hat{f}_z &= \sqrt{\frac{f}{2}} \left(\hat{\sigma}_{\uparrow \downarrow} + \hat{\sigma}_{\downarrow \uparrow}\right) + \sqrt{\frac{2f-1}{2}}  \left(\hat{\sigma}_{\downarrow T} + \hat{\sigma}_{T \downarrow }\right),
\end{align}
\end{subequations}
where we have defined the atomic population and coherence operators $\hat{\sigma}_{ba}=\ket{b}\bra{a}$.

\begin{figure}[htb]
\centering
  \includegraphics[width=.45\textwidth]{fig/FaradaySqueezingLevelStructure}
  \caption{The truncated qutrit subspace ground levels and related coherence transitions using a D2 line probe for cesium atoms in the optical pumping process.}\label{fig:spinsqueezinglevelstructure}
\end{figure}

Secondly, we assume the collective state is symmetric under exchange of spins. This approximation is valid when all atoms see the same probe intensity (i.e., they are sufficiently cold) and in the far-detuning regime when all sub- and super-radiant modes are equally excited. With this, we can limit our attention to the symmetric subspace and define, for example, the symmetric two-body covariances by
\begin{align}
\expect{\Delta\sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)}}_s \equiv \frac{1}{2}\left[\expect{\Delta\sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)}}+\expect{\Delta\sigma_{ba}^{(2)}\Delta\sigma_{dc}^{(1)}} \right] ,
\end{align}
where the superscripts, (1) and (2), label an arbitrary two atoms in the ensemble. Due to the exchange symmetry, $ \expect{\Delta\sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)}}_s=\expect{\Delta\sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)}}=\expect{\Delta\sigma_{ba}^{(2)}\Delta\sigma_{dc}^{(1)}} $, which reduces the number of $ n $-body moments required to simulate the spin dynamics of the ensemble.

Thirdly, we make the Gaussian approximation, valid for large atomic ensembles, so that the many-body state is fully characterized by one- and two-body correlations. Equivalently, the state is defined by the one- and two-body density operators, with matrix elements $\rho^{(1)}_{a, b} =\expect{\hat{\sigma}_{ba}}$, $\rho^{(1,2)}_{ac,bd}=\expect{\Delta \sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)} }_s$ in the symmetric subspace.   Optical pumping, acting locally, couple only $n$-body correlations to themselves, e.g.,
\begin{equation}\label{eq:dsigmabadc_op}
\left.d\expect{\Delta \sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)} }_s\right|_{op} = \expect{\mathcal{D}^\dagger[\Delta \sigma_{ba}^{(1)}]\Delta\sigma_{dc}^{(2)} }_sdt + \expect{\Delta \sigma_{ba}^{(1)} \mathcal{D}^\dagger[\Delta\sigma_{dc}^{(2)}] }_sdt .
\end{equation}
QND measurement generates higher order correlations according to
\begin{equation}\label{eq:dsigmaba_QND}
\left.d\expect{\hat{\sigma}_{ba}}\right|_{QND} =\frac{\kappa}{4}\expect{\mathcal{L}^\dagger\left[\hat{\sigma}_{ba} \right]}dt + \sqrt{\frac{\kappa}{4}}\expect{\mathcal{H}^\dagger\left[\hat{\sigma}_{ba} \right]}dW .
\end{equation}
We can truncate this hierarchy in the Gaussian approximation, setting third order cumulants to zero.  Thus, for example,
\begin{align}
\left.d\expect{\Delta \sigma_{ba}^{(1)} \Delta \sigma_{dc}^{(2)}}_s \right|_{QND} &= \left.d\expect{\hat{\sigma}_{ba}^{(1)} \hat{\sigma}_{dc}^{(2)}}_s \right|_{QND} - \left. \expect{\hat{\sigma}_{ba}} \right|_{QND} \left( \left.d\expect{\hat{\sigma}_{dc}} \right|_{QND}\right) \nonumber\\
&\quad - \left. \expect{\hat{\sigma}_{dc}} \right|_{QND} \left( \left.d\expect{\hat{\sigma}_{ba}} \right|_{QND}\right)
- \left.d\expect{\sigma_{ba}} \right|_{QND}\left.d\expect{\sigma_{dc}} \right|_{QND} \nonumber \\
&= -\kappa\expect{\Delta \sigma^{(1)}_{ba}  \Delta F_z }_s \expect{\Delta F_z \Delta \sigma_{dc}^{(2)} }_sdt,\label{eq:dsigmabadc_QND}
\end{align}
where we have employed the Ito calculus $dW^2 = dt$.
Note, by setting the third-order cumulants to zero,  the contribution of the $\mathcal{L}$ superoperator to dynamics of the two-body covariances vanishes,  $ \left.d\expect{\Delta \sigma_{ba}^{(1)} \Delta \sigma_{dc}^{(2)}}_s\right|_\mathcal{L} =\frac{\kappa}{4}\expect{\mathcal{L}^\dagger\left[\Delta\sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)} \right]}_sdt=0 $.  This indicates the events of no-photon detection under the Gaussian state approximation do not affect atom-atom correlations;  the measurement backaction and squeezing arise from the homodyne detection in the guided modes. 

Using all the approximations above, the collective spin dynamics can be efficiently calculated for the ensemble of qutrits (dimension  $d=3$) with $ d^2=9 $ equations for the one-body quantity, $ \expect{\hat{\sigma}_{ba}} $, and $ d^2(d^2+1)/2=45 $ equations for the two-body covariances, $ \expect{\Delta \sigma_{ba}^{(1)}\Delta\sigma_{dc}^{(2)} }_s $, in the symmetric subspace independent of the number of atoms.  With this formalism in hand, we can calculate the squeezing parameter, Eq.\eqref{eq:xi2Faraday}, as a function of time by finding time-dependent solutions for the one-body averages $\expect{\hat{f}_x}$ and  $\expect{\Delta f_z^2}$, and the two-body correlations $\expect{\Delta f_z^{(1)} \Delta f_z^{(2)}}$ and hence the collective quantities in the Wiener process.  The detailed approach to calculating the collective spin dynamics can be found in Appendix~\ref{Sec::opticalpumpinginrotatingframe}. 

In Fig.~\ref{fig:xi_rpfix_NA_t}, we plot the reciprocal of the spin squeezing parameter as a function of time and its peak as a function of atom number, $ N_A $, for both optical nanofiber and the square waveguide cases. By placing atoms $ 200 $nm away from the nanofiber surface ($ r'\!_\perp=1.8a $), $ 6.3 $dB of squeezing may be achievable for $ 2500 $ atoms. This is a substantial enhancement when compared to the birefringence protocol of spin squeezing  studied in~\cite{Qi2016}, which showed a peak squeezing  of $ 4.7 $dB with  similar experimental parameters. Using the square waveguide platform with the same number of atoms placed $150 $nm away from the surface of the square waveguide, our calculation yields $12.9$dB squeezing. As we have shown in Sec.~\ref{Sec::QNDandCooperativityTheory}, the square waveguide geometry enhances the anisotropic contrast of the two orthogonal guided modes and dramatically reduces the relative local intensity when the atoms are placed on the $ y $-axis. This results in a large cooperativity, higher peak spin squeezing, achieved in shorter time, and with a slower decay when using the square waveguide compared to the nanofiber. Figs.~\ref{fig:nanofiber_peakxi_NA_rp1d8a} and~\ref{fig:swg_peakxi_NA_rp1d} show how the peak squeezing scales with the number of trapped atoms. Using the nanofiber, as along as the atom number is  $>1400 $, the peak  squeezing  is larger than in the birefringence protocol with $ 2500 $ atoms using similar trapping parameters. For the square waveguide case, even if the number of trapped atoms is just above $ 800 $, $ >10 $dB of  squeezing is attainable. 

\begin{figure}[htb]
\centering
\pgfplotsset{xticklabel style={
        /pgf/number format/fixed,
        /pgf/number format/precision=2
},
scaled x ticks=false}
%\pgfplotsset{yticklabel style={text width=3em,align=right}}
 \begin{minipage}[h]{0.95\linewidth}
 %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
  \subfloat[h][$\xi^{-2}_{nanofiber}(t)$]{
    %\setlength\figureheight{0.3\textwidth}
    %\setlength\figurewidth{0.3\textwidth}
    \input{fig/nanofiber_xi_t_rp1d8a_NA2500.tex}
    \label{fig:nanofiber_xi_t_rp1d8a_NA2500}
    }
    \hfill
  \subfloat[h][$ \xi^{-2}_{peak,nanofiber}(N_A) $]{
      %\setlength\figureheight{0.3\textwidth}
      %\setlength\figurewidth{0.3\textwidth}
      \label{fig:nanofiber_peakxi_NA_rp1d8a}
      \input{fig/nanofiber_peakxi_NA_rp1d8a.tex}
      }
   \end{minipage}\vfill
   \begin{minipage}[h]{0.95\linewidth}
    %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
     \subfloat[h][$\xi^{-2}_{swg}(t)$]{
       %\setlength\figureheight{0.3\textwidth}
       %\setlength\figurewidth{0.3\textwidth}
       \input{fig/swg_xi_t_rp1d_NA2500.tex}
       \label{fig:swg_xi_t_rp1d_NA2500}
       }
       \hfill
     \subfloat[h][$ \xi^{-2}_{peak,swg}(N_A) $]{
         %\setlength\figureheight{0.3\textwidth}
         %\setlength\figurewidth{0.3\textwidth}
         \label{fig:swg_peakxi_NA_rp1d}
         \input{fig/swg_peakxi_NA_rp1d.tex}
         }
   \end{minipage}
   %\end{tabular}
\caption{The inverse spin-squeezing parameter. Eq.~\eqref{eq:xi2Faraday}, as a function of time and atom number. Subfigs.~\protect\subref*{fig:nanofiber_xi_t_rp1d8a_NA2500} and~\protect\subref*{fig:swg_xi_t_rp1d_NA2500} are plot $ \xi^{-2} $  as a function of time in units of the optical pumping rate $\gamma_{op}$, for the cylindrical nanofiber and square waveguide, respectively. Similarly, Subfigs.~\protect\subref*{fig:nanofiber_peakxi_NA_rp1d8a} and~\protect\subref*{fig:swg_peakxi_NA_rp1d} are the plots of $ \xi^{-2} $ as a function of $ N_A $ correspondingly. }\label{fig:xi_rpfix_NA_t}
\end{figure}

{\color{red} In the absence of any other noise, the cooperativity of atom-light coupling increases as the atoms are placed closer to the waveguide surface. Figs.\ref{fig:nanofiber_peakxi_rp_NA2500} and~\ref{fig:swg_peakxi_rp_NA2500} show the peak squeezing  as a function of $ r\!_\perp $ for both nanofiber and square waveguide geometries with $2500$ atoms. With the same setting, we also plot out the cooperativity, $ C_1 $, in the logarithm scale in Figs.\ref{fig:nanofiber_C1_y} and~\ref{fig:swg_C1_y} as a function of atom radial distance to the center of both waveguide geometries. We find that $C_1$ is proportional to $ e^{-\gamma r\!_\perp} $ and the peak squeezing scales as $ \sqrt{OD} $, where $ \gamma\approx 1.65/a $ for the nanofiber and $ \gamma\approx 2.14\times 2/w $ for the SWG with $2500$ atoms. The cooperativity of the SWG increases faster than the nanofiber as the atoms approach to the waveguide surface. 

In our simulations, we have used a probe light that is $\SI{4}{GHz}$ red-detuned to the $\ket{6S_{1/2}, f=4}\rightarrow \ \ket{6P_{3/2},f'=3}$ D2 line ($ \lambda\approx \SI{852}{nm} $) transition of cesium atoms. In principle, a probe that utilizes the D1 line transition may generate a slightly larger spin squeezing compared to the D2 line due to a stronger mode leakage to the trapped atoms at a given $ r'_\perp $ position. However, since the hyperfine splitting of the $ 6P_{1/2} $ excited state levels is much larger than that of $ 6P_{3/2} $ levels and not smaller enough than the $ f=4 $ and $ f=3 $ splitting of the $ 5S_{1/2} $ ground levels, a far-detuning probe to the D1 line transition might cause considerable spin dependent tensor light shifts or couplings to the ground or other excited levels of the cesium atoms. For the purpose of this paper, we decide to avoid the complexity of using D1 line signals and will provide more realistic designs and optimizations in our future publications. }

\newpage
\begin{figure}[htb]
\centering
%\pgfplotsset{yticklabel style={text width=3em,align=right}}
 \begin{minipage}[h]{0.95\linewidth}
 %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
  \subfloat[h][$\xi^{-2}_{peak,nanofiber}(r\!_\perp)$]{
    %\setlength\figureheight{0.3\textwidth}
    %\setlength\figurewidth{0.3\textwidth}
    \input{fig/nanofiber_peakxi_rp_NA2500.tex}
    \label{fig:nanofiber_peakxi_rp_NA2500}
    }
    \hfill
  \subfloat[h][$ C_1(r\!_\perp) $ of a nanofiber]{
      %\setlength\figureheight{0.3\textwidth}
      %\setlength\figurewidth{0.3\textwidth}
      \label{fig:nanofiber_C1_y}
      \input{fig/nanofiber_C1_y.tex}
      }
   \end{minipage}\vfill
   \begin{minipage}[h]{0.95\linewidth}
    %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
     \subfloat[h][$\xi^{-2}_{peak,swg}(r\!_\perp)$]{
       %\setlength\figureheight{0.3\textwidth}
       %\setlength\figurewidth{0.3\textwidth}
       \input{fig/swg_peakxi_rp_NA2500.tex}
       \label{fig:swg_peakxi_rp_NA2500}
       }
       \hfill
     \subfloat[h][$ C_1(r\!_\perp) $ of a SWG]{
         %\setlength\figureheight{0.3\textwidth}
         %\setlength\figurewidth{0.3\textwidth}
         \label{fig:swg_C1_y}
         \input{fig/swg_C1_y.tex}
         }
   \end{minipage}
   %\end{tabular}
\caption{Peak squeezing parameter (Subfigs.~\protect\subref*{fig:nanofiber_peakxi_rp_NA2500} and~\protect\subref*{fig:swg_peakxi_rp_NA2500}) and cooperativity at the optimal azimuthal tapping poistion (Subfigs.~\protect\subref*{fig:nanofiber_C1_y} and~\protect\subref*{fig:swg_C1_y}) as a function of the radial distance to the waveguide axis. \protect\subref*{fig:nanofiber_peakxi_rp_NA2500} and~\protect\subref*{fig:nanofiber_C1_y} are for the nanofiber case; Subfigs.~\protect\subref*{fig:swg_peakxi_rp_NA2500} and~\protect\subref*{fig:swg_C1_y} are for the square waveguide case. }\label{fig:peakxi_rp_NA}
\end{figure}


%=================== Comparison and generalization =====================%
\subsection{Spin squeezing with rectangular waveguide and general criteria in optimizing spin squeezing effect} \label{Sec::Waveguide}

Generalized OD/$N_A$ using Green function method (or may be just with two orthogonal mode bases).

Theoretical upper bounds for the Birefringence and Faraday protocols using a waveguide:

1. The maximum achievable OD$ /N_A $ can be defined when the two orthogonal $ H $- and $ V $-modes are purely linearly polarized (transverse modes).

2. The relationship between the upper bound OD$ /N_A $ and spin-squeezing parameter due to the waveguide effect and effective mode area at the atom position. 
General rules to find the optimal choice of atom positions for the Birefringence and Faraday spin squeezing protocols.


As one example of near-linearly polarized modes, we can consider using a dielectric waveguide with a square intersection. 
Coupling strength, decoherence, and spin squeezing analysis based on the modes of the squared waveguides to show improvement compared with the nanofiber geometry...

Relation of OD$ /N_A $ to cooling efficiency and fictitious magnetic field applications.


\section{A two-color spin squeezing protocol to cancel the tensor light shift}
As illustrated in Enrique Montano's PhD dissertation work, the tensor light shift due to the external field can be canceled in a free-space spin squeezing setup, we will generalize this idea to the nanophotonic waveguide case.

From Eq.~\eqref{eq:Heff_Faraday_C02}, the state-dependent tensor light shift (the term proportional to $ \hat{f}_x^2 $) is proportional to 
\begin{align}
\delta E_T &= \sum_{f'} \frac{\Gamma_{f'}^x\Omega^2}{\Delta_{ff'}^2+(\Gamma_{f'}^x)^2/4}C_{ff'}^{(2)}\\
&\approx \sum_{f'}\frac{\Gamma_{f'}^x\Omega^2}{\Delta_{ff'}^2}C_{ff'}^{(2)} = \sum_{f'}\frac{\sigma_0}{\Ain}\left(\frac{\Gamma_{f'}^x}{\Delta_{ff'}} \right)^2 \dot{N}_L C_{ff'}^{(2)}
\end{align}



%====== SECTION: Summary and outlook ======%
\section{Summary and Outlook} \label{Sec::Conclusion}


ACKNOWLEDGMENTS
We thank the UNM Center for Advanced Research Computing for computational resources used in this work. We also thank Perry Rice and Ezad Shojaee for stimulated discussions. 
This work was supported by the NSF, under grants PHY-1212445, xxxxx.

\bibliography{refs/Archive}


%=========== APPENDIX ===========%
\begin{appendix}

%===================APPENDIX: Hamiltonian =====================%
\section{Faraday interaction Hamiltonian} \label{Appendix::FaradayInteractionHamiltonian}
In the Faraday interaction spin squeezing protocol, we define the fiducial, coupled and transfer states by $ \ket{\uparrow}=\ket{f=4,f_x=4} $, $ \ket{\downarrow}=\ket{f=4,f_x=3} $ and $ \ket{T}=\ket{f=4,f_x=2} $ respectively in the $ x $-basis. 
A set of spin operators projected onto the truncated qutrit subspace spanned by these three basis states can be defined by
\begin{align}
\hat{f_x} &= -f \ket{\uparrow}\bra{\uparrow} -(f-1)\ket{\downarrow}\bra{\downarrow}-(f-2)\ket{T}\bra{T},\\
\hat{f_y} &=i\left[\sqrt{\frac{f}{2}}\left(\ket{\downarrow}\bra{\uparrow}-\ket{\uparrow}\bra{\downarrow}\right) +\sqrt{\frac{2f-1}{2}}\left(\ket{T}\bra{\downarrow}-\ket{\downarrow}\bra{T} \right) \right] ,\\
\hat{f_z} &= \sqrt{\frac{f}{2}}\left(\ket{\downarrow}\bra{\uparrow}+\ket{\uparrow}\bra{\downarrow}\right) +\sqrt{\frac{2f-1}{2}}\left(\ket{T}\bra{\downarrow}+\ket{\downarrow}\bra{T} \right).
\end{align}

We also define a set of Stokes operators as below,
\begin{align}
\hat{S}_0 &= \smallfrac{1}{2}\big[ \hat{a}^\dag_H(t) \hat{a}_H(t)+\hat{a}^\dag_V(t) \hat{a}_V(t) \big],\\
\hat{S}_1 &= \smallfrac{1}{2}\big[ \hat{a}^\dag_H(t) \hat{a}_H(t)-\hat{a}^\dag_V(t) \hat{a}_V(t) \big],\\
\hat{S}_2 &= \smallfrac{1}{2}\big[ \hat{a}^\dag_H(t) \hat{a}_V(t)+\hat{a}^\dag_V(t) \hat{a}_H(t) \big],\\
\hat{S}_3 &= \smallfrac{1}{2i}\big[ \hat{a}^\dag_H(t) \hat{a}_V(t) -\hat{a}^\dag_V(t) \hat{a}_H(t) \big],
\end{align}
where the photon annihilation operators of the horizontally and vertically polarized modes denoted with subscription $ H $ and $ V $ respectively are related to the left- and right-circularly polarized modes denoted with subscription $ L $ and $ R $ by $ \hat{a}_{H}=(\hat{a}_L+\hat{a}_R) /\sqrt{2}$ and $ \hat{a}_{V}=i(\hat{a}_R-\hat{a}_L)/\sqrt{2} $.
The local quasimonochromatic electric field operator at $\br= (r\!_\perp,\phi,z) $ with a set of degenerate guided modes at frequency $ \omega_0 $ and propagating in the group velocity of $ v_g $ can be given by
\begin{align}\label{eq:Ebp}
\hat{\mathbf{E}}^{(+)}(r\!_\perp,\phi,z;t) &= \sum_{b,p} \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \mathbf{u}_{b,p}(r\!_\perp,\phi) \hat{a}_{b,p}(z,t)  e^{i (b\beta_0 z- \omega_0 t)},
\end{align}
where $ b $ is the sign of the propagation direction and $ p $ is the polarization label of the orthonormal modes.
In our case, we only have forward propagating modes for the QND measurement protocols, and hence we will only have $ b=+ $.

\comment{Todo: Rewrite the effective Hamiltonian using the collective operators defined above.}

%===================APPENDIX: Circular VS linear polarization mode bases =====================%
\section{Circular V.S. linear polarization bases}\label{Appendix:LRbases}
\comment{This content is mainly typed from my notes on Faraday Interaction in The Nanofiber Geometry. More can be done here, including a generalization of polarization transformation using $ \mathbf{SU}(2) $ group generators.}

In this section, we will derive the equations of Stokes vector operators and spin-polarization interaction Hamiltonian due to polarization basis transformations.
A general basis transformation theory will be derived in the context of linear-polarization basis transformations, and will then be applied to the linear $ D/\bar{D} $- and circular $ R/L $-bases cases.

\subsection{Spin-polarization Hamiltonian and Stokes vector operators in a linear basis}
In general, we can define an arbitrary linear polarization basis by
\begin{align}
\left(\!\begin{array}{c}
\mathbf{e}_n \\ \mathbf{e}_{\bar{n}}
\end{array}\!\right) &= 
\left(\!\!\begin{array}{cc}
\cos\theta & \sin\theta \\
- \sin\theta & \cos\theta
\end{array}\!\!\right)\bullet
\left(\!\begin{array}{c}
\mathbf{e}_H \\ \mathbf{e}_V
\end{array}\!\right)
=\mathbf{R}(\theta)\bullet \left(\!\begin{array}{c}
\mathbf{e}_H \\ \mathbf{e}_V
\end{array}\!\right),
\end{align}
or the inversed relationship
\begin{align}
\left(\!\begin{array}{c}\mathbf{e}_H \\ \mathbf{e}_V\end{array}\!\right)&= \mathbf{R}^{-1}(\theta)\bullet\left(\!\begin{array}{c}\mathbf{e}_n \\ \mathbf{e}_{\bar{n}}\end{array}\!\right),
\end{align}
where $ \theta $ is the angle of the $ \mathbf{e}_n $ basis rotated from the $ H $-direction around $ z $-axis, and $ \mathbf{e}_{\bar{n}} $ is the basis vector $ 90^\circ $ from the $ \mathbf{e}_n $ direction; $ \mathbf{R}(\theta)=\mathbf{R}_z(\theta) $ is the Euler rotation matrix about the $ z $-axis by $ \theta $ in the real-number $ \mathbf{SO}(3) $ rotation group, which has the property that $ \mathbf{R}^{-1}(\theta)=\mathbf{R}^T(\theta)=\mathbf{R}(-\theta) $. 
More generally, the basis transformation matrix is an unitary matrix determined by two parameters (two degrees of freedom)--$ \theta $ and $ \phi $--corresponding to the rotating angles around one axis and an relative phase between the base components, which is in the $ \mathbf{SU}(2) $ group.
%We denote the general case with $ \mathbf{R}=\mathbf{R}(\theta,\phi) $, or in the form of two-step rotations around $ i $-axis and then around $ j $-axis by $ \mathbf{R}=\mathbf{R}_j(\theta)\mathbf{R}_i(\phi) $, always satisfying $ \mathbf{R}^{-1}=\mathbf{R}^\dagger $ for either rotations.

Not to be confused, we have also defined an operator space spanned by operator vectors, like $ \left(\!\begin{array}{cc}\mathbf{e}_n,&\mathbf{e}_{\bar{n}}\end{array}\! \right) $, which has vectors, tensors or operators as the elements.
We have also defined the bullet operator ($ \bullet $) in the operator vector space isomorphically the same as the dot ($ \cdot $) product or matrix product in the conventional vector space while the sign of $ \cdot $ can usually be ignored and we will denote complex conjugates explicitly if needed. 
$ \mathbf{R}(\theta) $ and its transformations is a tensor defined in the operator vector space as well.
When a conventional vector or tensor multiplies with an operator vector or tensor, we will use $ \cdot $ between them and the conventional vector or tensor will be formally treated as a scalar to be $ \cdot $ multiplied with the elements of the operator vector or tensor. 
Two operator vectors in a $ \bullet $ multiplication form a mutual covariant relationship in the operator space.
 
With the coordinate basis rotated passively, both the mode components and the field annihilation operators should be rotated actively by $ -\theta $ to be transformation-equivalent.
Written in the matrix form in the operator space, 
\begin{align}
\left(\!\begin{array}{c}\mathbf{u}_n \\ \mathbf{u}_{\bar{n}}\end{array}\!\right) &= \mathbf{R}^{-1}(\theta)\bullet\left(\!\begin{array}{c}\mathbf{u}_H \\ \mathbf{u}_V\end{array}\!\right)\\
\hat{\mathbf{a}}_{n,\bar{n}} &=\mathbf{R}^{-1}(\theta)\bullet \hat{\mathbf{a}}_{H,V},
\end{align}
where $ \hat{\mathbf{a}}_{n,\bar{n}}=[\hat{a}_n;\hat{a}_{\bar{n}}] $ and $ \hat{\mathbf{a}}_{H,V}=[\hat{a}_H;\hat{a}_V] $ are the annihilation operator vectors using the $ \{n,\bar{n} \} $ and $ \{H,V \} $ bases, respectively.
In our notation, we use $ [\cdot ;\cdots] $ notation to indicate $ 1\times n $ vectors as general matrices.

Using the definition in the $ \{ H,V\} $-basis of the Stokes operators (Eq.\eqref{eq:SaHVaRL}) and the annihilation operator basis transformation relationships above, one can rewrite the Stokes operators in the $ \{\mathbf{e}_n, \mathbf{e}_{\bar{n}}\} $ basis by
\begin{subequations}\label{eq:Snnbar}
\begin{align}
\hat{S}_0 &= \frac{1}{2}\left[\left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)\bullet\mathbf{R}_{[1,:]}^\dagger(\theta)\bullet\mathbf{R}_{[1,:]}(\theta)\bullet
\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right)
+ \left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)\bullet\mathbf{R}^\dagger_{[2,:]}(\theta)\bullet\mathbf{R}_{[2,:]}(\theta)\bullet
\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \right]\nn\\
&= \frac{1}{2}\left\{\left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)\bullet
\left[\left(\!\begin{array}{cc}\cos^2\theta,& \frac{1}{2}\sin 2\theta \\ \frac{1}{2}\sin 2\theta, & \sin^2\theta\end{array} \!\right)
+ \left(\!\begin{array}{cc}\sin^2\theta,& -\frac{1}{2}\sin 2\theta \\ -\frac{1}{2}\sin 2\theta, & \cos^2\theta\end{array} \!\right)\right]\bullet
\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right)\right\}\nn\\
&=\frac{1}{2} \left[\hat{a}_n^\dagger\hat{a}_n+\hat{a}_{\bar{n}}^\dagger\hat{a}_{\bar{n}} \right]\\
\hat{S}_1 &= \frac{1}{2}\left\{\left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)
\bullet\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\bullet\mathbf{R}_{[1,:]}(\theta) - \mathbf{R}_{[2,:]}^\dagger(\theta)\bullet\mathbf{R}_{[2,:]}(\theta) \right]
\bullet\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \right\}\nn\\
&= \frac{1}{2} \left[\cos 2\theta \hat{a}_n^\dagger\hat{a}_n + \sin 2\theta \hat{a}_n^\dagger\hat{a}_{\bar{n}} + \sin 2\theta \hat{a}_{\bar{n}}^\dagger\hat{a}_n -\cos 2\theta \hat{a}_{\bar{n}}^\dagger\hat{a}_{\bar{n}} \right]\\
\hat{S}_2 &= \frac{1}{2}\left\{\left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)
\bullet\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\bullet\mathbf{R}_{[2,:]}(\theta) + \mathbf{R}_{[2,:]}^\dagger(\theta)\bullet\mathbf{R}_{[1,:]}(\theta) \right]
\bullet\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \right\}\nn\\
&= \frac{1}{2} \left[-\sin 2\theta \hat{a}_n^\dagger\hat{a}_n + \cos 2\theta \hat{a}_n^\dagger\hat{a}_{\bar{n}} + \cos 2\theta \hat{a}_{\bar{n}}^\dagger\hat{a}_n +\sin 2\theta \hat{a}_{\bar{n}}^\dagger\hat{a}_{\bar{n}} \right]\\
\hat{S}_3 &= \frac{1}{2i}\left\{\left(\!\begin{array}{cc}\hat{a}_n^\dagger,& \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)
\bullet\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\bullet\mathbf{R}_{[2,:]}(\theta) - \mathbf{R}_{[2,:]}^\dagger(\theta)\bullet\mathbf{R}_{[1,:]}(\theta) \right]
\bullet\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \right\}\nn\\
&= \frac{1}{2i} \left[\hat{a}_n^\dagger\hat{a}_{\bar{n}} - \hat{a}_{\bar{n}}^\dagger\hat{a}_n  \right].
\end{align}
\end{subequations}
In deriving the equations above, we have denoted $ \mathbf{R}_{[i,:]}(\theta) $ as the $ i $-th row of $ \mathbf{R}(\theta) $ and $ \mathbf{R}_{[i,:]}^\dagger(\theta) $ as the conjugate transpose of the $ i $-th row of $ \mathbf{R}(\theta) $.
As a shorthand, these relationships of Stokes operators can be expressed as an operator transformation, $ \hat{\mathbf{S}}=\mathbf{M}\bullet\hat{\mathbf{A}}_{n,\bar{n}} $, where $ \hat{\mathbf{S}}=[\hat{S}_0;\hat{S}_1;\hat{S}_2;\hat{S}_3] $ and $ \hat{\mathbf{A}}_{n,\bar{n}}=[\hat{a}_n^\dagger\hat{a}_n;\hat{a}_n^\dagger\hat{a}_{\bar{n}};\hat{a}_{\bar{n}}^\dagger\hat{a}_n;\hat{a}_{\bar{n}}^\dagger\hat{a}_{\bar{n}}] $ are the operator vectors and $ \mathbf{M} $ is the transformation matrix defined by the transformation coefficients in Eqs.\eqref{eq:Snnbar}.
One can prove that the inversed transformation matrix $ \mathbf{M}^{-1}=2\mathbf{M}^\dagger $, and $ \mathbf{M} $ can be derived in the following form, in general,
\begin{align}
\mathbf{M} &=\frac{1}{2}\left(\!\begin{array}{c}
\mathrm{vec}_r\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\mathbf{R}_{[1,:]}(\theta) + \mathbf{R}_{[2,:]}^\dagger(\theta)\mathbf{R}_{[2,:]}(\theta) \right]\\
\mathrm{vec}_r\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\mathbf{R}_{[1,:]}(\theta) - \mathbf{R}_{[2,:]}^\dagger(\theta)\mathbf{R}_{[2,:]}(\theta) \right]\\
\mathrm{vec}_r\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\mathbf{R}_{[2,:]}(\theta) + \mathbf{R}_{[2,:]}^\dagger(\theta)\mathbf{R}_{[1,:]}(\theta) \right]\\
-i\mathrm{vec}_r\left[\mathbf{R}_{[1,:]}^\dagger(\theta)\mathbf{R}_{[2,:]}(\theta) - \mathbf{R}_{[2,:]}^\dagger(\theta)\mathbf{R}_{[1,:]}(\theta) \right]
 \end{array} \!\right)\\
&= \frac{1}{2}\left(\!\begin{array}{cccc} 1,&0,&0,& 1\\
\cos 2\theta,&\sin 2\theta, & \sin 2\theta, & -\cos 2\theta\\
-\sin 2\theta, & \cos 2\theta, & \cos 2\theta, & \sin 2\theta\\
0,&-i,& i,& 0\end{array}\!\right),
\end{align}
where $ \mathrm{vec}_r[\cdot] $ means the vectorization of a matrix by concatenating its rows. 


By using the vector space representation, the E-field operator can be written in the $ \{n,\bar{n}\} $ basis by 
\begin{align}
\hat{\mathbf{E}}^{(+)}(\br;t) &= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left[\mathbf{u}_H(r\!_\perp,\phi) \hat{a}_H(z,t) + \mathbf{u}_V(r\!_\perp,\phi) \hat{a}_V(z,t)\right]  e^{i (\beta_0 z- \omega_0 t)}\\
&= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left(\!\begin{array}{cc}\mathbf{u}_H, & \mathbf{u}_V\end{array}\!\right)
\bullet \left(\!\begin{array}{c}\hat{a}_H\\ \hat{a}_V\end{array} \!\right)  e^{i (\beta_0 z- \omega_0 t)}\\
&= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left(\!\begin{array}{cc}\mathbf{u}_H, & \mathbf{u}_V\end{array}\!\right)\bullet\mathbf{R}(\theta)
\bullet \mathbf{R}^{-1}(\theta)\bullet\left(\!\begin{array}{c}\hat{a}_H\\ \hat{a}_V\end{array} \!\right)  e^{i (\beta_0 z- \omega_0 t)}\\
&= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left(\!\begin{array}{cc}\mathbf{u}_n, & \mathbf{u}_{\bar{n}}\end{array}\!\right)
\bullet \left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right)  e^{i (\beta_0 z- \omega_0 t)}\\
&= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left[\mathbf{u}_n(r\!_\perp,\phi) \hat{a}_n(z,t) + \mathbf{u}_{\bar{n}}(r\!_\perp,\phi) \hat{a}_{\bar{n}}(z,t)\right]  e^{i (\beta_0 z- \omega_0 t)}.\label{eq:Efieldop_nnbar}
\end{align}
As expected, the field operator preserves the form in the new basis as defined in Eq.\eqref{eq:Ebp}.

Therefore, the effective atom-light interaction Hamiltonian can be written as
\begin{align}
\hat{h}_\eff &= -\hat{\mathbf{E}}^{(-)}(\br')\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot\hat{\mathbf{E}}^{(+)}(\br')\nn\\
&=-\frac{2 \pi \hbar \omega_0}{ v_g}
\left(\!\begin{array}{cc}\hat{a}_n^\dagger, & \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)\bullet \left(\!\begin{array}{c}\mathbf{u}_n^*\\ \mathbf{u}_{\bar{n}}^*\end{array}\!\right)
\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot
\left(\!\begin{array}{cc}\mathbf{u}_n, & \mathbf{u}_{\bar{n}}\end{array}\!\right)
\bullet \left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \nn\\
&= \hbar \left(\!\begin{array}{cc}\hat{a}_n^\dagger, & \hat{a}_{\bar{n}}^\dagger\end{array} \!\right)\bullet
\left(\!\begin{array}{cc} \hat{\chi}_{nn},&\hat{\chi}_{n\bar{n}}\\
\hat{\chi}_{\bar{n}n},&\hat{\chi}_{\bar{n}\bar{n}}\end{array} \!\right)
\bullet\left(\!\begin{array}{c}\hat{a}_n\\ \hat{a}_{\bar{n}}\end{array} \!\right) \nn\\
&= \hbar \hat{\boldsymbol{\chi}}_{n,\bar{n}}^\dagger\bullet \hat{\mathbf{A}}_{n,\bar{n}}= \hbar \hat{\boldsymbol{\chi}}_{n,\bar{n}}^\dagger\bullet 2\mathbf{M}^\dagger\bullet \hat{\mathbf{S}}\\
&= \hbar\left\{(\hat{\chi}_{nn}+ \hat{\chi}_{\bar{n}\bar{n}})\hat{S}_0 \right.\nn\\
&\quad+ [\cos 2\theta \hat{\chi}_{nn}+ \sin 2\theta(\hat{\chi}_{n\bar{n}}+\hat{\chi}_{\bar{n}n}) - \cos 2\theta\hat{\chi}_{\bar{n}\bar{n}}]\hat{S}_1 \nn\\
&\quad+ [-\sin 2\theta \hat{\chi}_{nn}+ \cos 2\theta(\hat{\chi}_{n\bar{n}}+\hat{\chi}_{\bar{n}n}) + \sin 2\theta \hat{\chi}_{\bar{n}\bar{n}}]\hat{S}_2 \nn\\
&\quad+\left. i \left(\hat{\chi}_{n\bar{n}}-\hat{\chi}_{\bar{n}n}\right)\hat{S}_3 \right\}.\label{eq:heff_nnbarChiS}
\end{align}
In deriving the expression above, we have defined $\hat{\boldsymbol{\chi}}_{n,\bar{n}}=[\hat{\chi}_{nn};\hat{\chi}_{n\bar{n}};\hat{\chi}_{\bar{n}n};\hat{\chi}_{\bar{n}\bar{n}}]$ as the coupling operator vector.
The coupling coefficient elements of $ 2\mathbf{M}^\dagger $, $ 2M_{ij} $, correspond to the coupling coefficients in Eq.\eqref{eq:heff_nnbarChiS} between the spin operators in $ \hat{\boldsymbol{\chi}}_{n,\bar{n}} $ and the polarization operators in $ \hat{\mathbf{S}} $, and each column of $ 2\mathbf{M}^\dagger $ corresponds to the coefficients in the corresponding line of the $ \hat{S}_i $ coupling term in Eq.\eqref{eq:heff_nnbarChiS}.


If we set $ \theta=\pi/4 $, we will be in the $ \{\mathbf{e}_D,\mathbf{e}_{\bar{D}} \} $ basis.
\begin{subequations}
\begin{align}
\hat{a}_D &= \frac{1}{\sqrt{2}}(\hat{a}_H-\hat{a}_V)\\
\hat{a}_{\bar{D}} &= \frac{1}{\sqrt{2}}(\hat{a}_H+\hat{a}_V),
\end{align}
\end{subequations}
or,
\begin{subequations}
\begin{align}
\hat{a}_H &= \frac{1}{\sqrt{2}}(\hat{a}_D+\hat{a}_{\bar{D}})\\
\hat{a}_V &= \frac{1}{\sqrt{2}}(\hat{a}_{\bar{D}}-\hat{a}_D).
\end{align}
\end{subequations}
The Stokes operators becomes
\begin{subequations}
\begin{align}
\hat{S}_0 &= \frac{1}{2} \left[\hat{a}_D^\dagger\hat{a}_D+\hat{a}_{\bar{D}}^\dagger\hat{a}_{\bar{D}} \right]\\
\hat{S}_1 &= \frac{1}{2} \left[ \hat{a}_D^\dagger\hat{a}_{\bar{D}} +  \hat{a}_{\bar{D}}^\dagger\hat{a}_D \right]\\
\hat{S}_2 &= \frac{1}{2} \left[ \hat{a}_{\bar{D}}^\dagger\hat{a}_{\bar{D}}-\hat{a}_D^\dagger\hat{a}_D \right]\\
\hat{S}_3 &= \frac{1}{2i} \left[\hat{a}_D^\dagger\hat{a}_{\bar{D}}-\hat{a}_{\bar{D}}^\dagger\hat{a}_D \right].
\end{align}
\end{subequations}
The field operator becomes 
\begin{align}
\hat{\mathbf{E}}^{(+)}(\br;t) &=\sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left[\mathbf{u}_D(r\!_\perp,\phi) \hat{a}_D(z,t)+\mathbf{u}_{\bar{D}}(r\!_\perp,\phi) \hat{a}_{\bar{D}}(z,t)  \right]e^{i (\beta_0 z- \omega_0 t)}.
\end{align}
The effective Hamiltonian can be simplified from Eq.\eqref{eq:heff_nnbarChiS} to 
\begin{align}
\hat{h}_\eff &= \hbar[(\hat{\chi}_{DD}+\hat{\chi}_{\bar{D}\bar{D}})\hat{S}_0 +(\hat{\chi}_{D\bar{D}}+\hat{\chi}_{\bar{D}D} )\hat{S}_1\nn\\
&\quad + (\hat{\chi}_{\bar{D}\bar{D}}-\hat{\chi}_{DD})\hat{S}_2 +i(\hat{\chi}_{D\bar{D}}-\hat{\chi}_{\bar{D}D} )\hat{S}_3].
\end{align}
The Hamiltonian expression above should satisfy the cyclical transformation from the $ H $- and $ V $-basis expression.

\subsection{Spin-polarization Hamiltonian and Stokes vector operators in a circular basis}
We define a set of polarization vector transformation relationships by 
\begin{subequations}
\begin{align}
\hat{a}_H &= \frac{1}{\sqrt{2}}(\hat{a}_R+\hat{a}_L )\\
\hat{a}_V &= \frac{i}{\sqrt{2}}(\hat{a}_R-\hat{a}_L ),
\end{align}
\end{subequations}
or the inverse
\begin{subequations}
\begin{align}
\hat{a}_R &= \frac{1}{\sqrt{2}}(\hat{a}_H-i\hat{a}_V )\\
\hat{a}_L &= \frac{1}{\sqrt{2}}(\hat{a}_H+i\hat{a}_V ),
\end{align}
\end{subequations}
where $ R $($ L $) indicates the right(left)-circularly polarized mode.
The Stokes operators can then be defined in both linear ($ H $ and $ V $) and circular ($ L $ and $ R $) polarization bases by
\begin{subequations}\label{eq:SaHVaRL}
\begin{align}
\hat{S}_0 &= \frac{1}{2} \left[\hat{a}_H^\dagger\hat{a}_H+\hat{a}_V^\dagger\hat{a}_V \right] = \frac{1}{2} \left[\hat{a}_R^\dagger\hat{a}_R+\hat{a}_L^\dagger\hat{a}_L \right]\\
\hat{S}_1 &= \frac{1}{2} \left[\hat{a}_H^\dagger\hat{a}_H-\hat{a}_V^\dagger\hat{a}_V \right] = \frac{1}{2} \left[\hat{a}_R^\dagger\hat{a}_L+\hat{a}_L^\dagger\hat{a}_R \right]\\
\hat{S}_2 &= \frac{1}{2} \left[\hat{a}_H^\dagger\hat{a}_V+\hat{a}_V^\dagger\hat{a}_H \right] = \frac{i}{2} \left[\hat{a}_L^\dagger\hat{a}_R-\hat{a}_R^\dagger\hat{a}_L \right]\\
\hat{S}_3 &= \frac{1}{2i} \left[\hat{a}_H^\dagger\hat{a}_V-\hat{a}_V^\dagger\hat{a}_H \right] = \frac{1}{2} \left[\hat{a}_R^\dagger\hat{a}_R-\hat{a}_L^\dagger\hat{a}_L \right].
\end{align}
\end{subequations}
The inversed transformations can be easily derived by inverting the transformation coefficient matrices. 

Based on Eq.\eqref{eq:Ebp}, when the input probe can be decomposed into degenerate orthonormal guided modes in the $H/V $ or $ R/L $ bases, the E-field operator can be written in the quasilinear and quasicircular mode bases by 
\begin{align}
\hat{\mathbf{E}}^{(+)}(r\!_\perp,\phi,z;t) &= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left[\mathbf{u}_H(r\!_\perp,\phi) \hat{a}_H(z,t) + \mathbf{u}_V(r\!_\perp,\phi) \hat{a}_V(z,t)\right]  e^{i (\beta_0 z- \omega_0 t)}\\
&= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g} } \left[\mathbf{u}_R(r\!_\perp,\phi) \hat{a}_R(z,t) + \mathbf{u}_L(r\!_\perp,\phi) \hat{a}_L(z,t)\right]  e^{i (\beta_0 z- \omega_0 t)}.
\end{align}
Therefore, the effective atom-light interaction Hamiltonian can be given in those bases by 
\begin{align}
\hat{h}_\eff &= -\hat{\mathbf{E}}^{(-)}(\br')\cdot\hat{\tensor{\mathbf{\alpha}}}\cdot\hat{\mathbf{E}}^{(+)}(\br')\nn\\
&= \hbar\left[(\hat{\chi}_{HH}+\hat{\chi}_{VV})\hat{S}_0 + (\hat{\chi}_{HH}-\hat{\chi}_{VV})\hat{S}_1 + (\hat{\chi}_{HV}+\hat{\chi}_{VH})\hat{S}_2 + i(\hat{\chi}_{HV}-\hat{\chi}_{VH})\hat{S}_3 \right]\\
&= \hbar\left[(\hat{\chi}_{RR}+\hat{\chi}_{LL})\hat{S}_0 + (\hat{\chi}_{RL}+\hat{\chi}_{LR})\hat{S}_1 + i(\hat{\chi}_{RL}-\hat{\chi}_{LR})\hat{S}_2 + (\hat{\chi}_{RR}-\hat{\chi}_{LL})\hat{S}_3 \right]\\
&=\hbar\sum_{i=0}^3 \hat{\chi}_{i}\hat{S}_i\\
&=\hbar\sum_{i,j=0} \chi_{ij}\hat{f}_i\hat{S}_j,
\end{align}
with $\hat{\chi}_{pp'} $ defined in Eq.\eqref{eq:chippp}.


%===================APPENDIX: Rectangular waveguides and mode walkoff ==================== %
\section{Rectangular waveguides and mode walkoff}
Based on Eq.\eqref{eq:Ebp}, if there are two non-degenerate guided modes used for the QND measurement, we will have
\begin{align}\label{eq:EdifferentHV}
\hat{\mathbf{E}}^{(+)}(r\!_\perp,\phi,z;t) &= \sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g^H} } \mathbf{u}_H(r\!_\perp,\phi) \hat{a}_H(z,t)  e^{i (\beta_0^H z- \omega_0 t)}\nn\\
&\quad +\sqrt{ \frac{2 \pi \hbar \omega_0}{ v_g^V} } \mathbf{u}_V(r\!_\perp,\phi) \hat{a}_V(z,t)  e^{i (\beta_0^V z- \omega_0 t)},
\end{align}
where $ v_g^{H/V} $ and $ \beta_0^{H/V} $ are the group velocities and propagation constants of the $ H $- and $ V $-modes, respectively. 
The offset of $ v_g^{H/V} $ and $ \beta_0^{H/V} $ for the two guided modes could lead to a walkoff effect of the modes and cause an intrinsic phase shift.

%===================APPENDIX: Photon scattering and optical pumping rates =====================%
\section{Photon scattering and optical pumping rates} \label{Appendix::Rates}

In this Appendix, we give the explicit expressions for the photon scattering rates following the formalism outlined in~\cite{Deutsch2010a} and ~\cite{Qi2016} while also considering the modification of decay rates of atoms due to the presence of the nanophotonic structures. 

We define the Green's function tensor in presence of a dielectric medium satisfying the equation~\cite{Qi2016},
	\begin{align} \label{Eq::GreensDiffEq}
		\left[ -\nabla\times\nabla\times + n^2(\mbf{r}) k_0^2 \right] \tensor{\mathbf{G}}(\br, \br';\omega_0) &= -4\pi 
k_0^2 \delta^{(3)}(\mathbf{r}-\mathbf{r}') \unittens,
	\end{align}
where $\unittens$ is the unit tensor.
Under the Born approximation, the field at $\br$ due to an input field $ \mathbf{E}_{\inp}(\br) $ and a dipole source placed at $\br'$ can then be written as
\begin{align}
		\mathbf{E}_{\out}(\br) 
		&\approx \mathbf{E}_{\inp}(\br)+ \tensor{\mathbf{G}}^{(+)}(\br , \br'; \omega_0) \cdot 
\tensor{\boldsymbol{\alpha}}\cdot \mathbf{E}_{\inp}(\br'), \label{Eq::ScatteredField}
\end{align}
where $ \tensor{\boldsymbol{\alpha}} $ is the polarizability of the dipole source.
We can specify the dipole source as an atom with arbitrary excited state level $ e $ and ground level $ g $ in the hyperfine structure manifold. 
The atomic decay rate of the excited level $ e $ will be modified due to the presence of the nanophotonic structure and can be given by
\begin{align}
\Gamma_e &= \frac{2}{\hbar}\sum_g \bra{g}\hat{\mathbf{d}}\ket{e} \cdot \mathrm{Im}\left[\tensor{\mathbf{G}}(\br',\br';\omega_{eg}) \right]\cdot \bra{e}\hat{\mathbf{d}}\ket{g}\label{eq:Gamma_e}
\end{align}
with the total Green's function tensor as a combination of the contributions from the guided modes and radiative modes, $ \tensor{\mathbf{G}}(\br',\br';\omega_{eg})=\tensor{\mathbf{G}}_{gyd}(\br',\br';\omega_{eg})+\tensor{\mathbf{G}}_{rad}(\br',\br';\omega_{eg}) $.
For simplicity, we will hide the frequency-dependence part from the Green's function tensors below.
Correspondingly, $ \Gamma_e=\Gamma_{e,gyd}+\Gamma_{e,rad} $ can also be decomposed into the guided mode and radiative mode contribution parts. 
In the free space, the local Green's function tensor has a divergent real part and the imaginary part in the CGS units can be given by 
$\mathrm{Im}\left[\tensor{\mathbf{G}}_0\right]=G_0\tensor{\mathbbm{1}}$ with $G_0(\mathbf{r}',\mathbf{r}';\omega_{eg})=\frac{2}{3}k_0^3$.
In presence of a waveguide, the Green's function tensor due to the radiative modes can be decomposed into the free-space or homogeneous radiation contribution and the waveguide-induced radiation contribution, that is $ \mathrm{Im}\left[\tensor{\mathbf{G}}_{rad}(\br',\br')\right]=\tensor{\mathbf{G}}_0(\br',\br')+\tensor{\mathbf{G}}_{ind,rad}(\br',\br';\omega_{eg}) $.
The waveguide-induced Green's function tensor elements due to the radiative modes can be calculated by
$$G_{ind,rad}^{ij}(\mathbf{r}',\mathbf{r}')=\frac{2k_0^2}{3\pi^2}\int_{-k_0}^{k_0} d\beta E_j^i(\mathbf{r}')=\frac{4k_0^2}{3\pi^2}\int_{0}^{k_0} d\beta E_j^i(\mathbf{r}')$$
once $ E_j^i(\mathbf{r}') $--the $i$-th electric field component measured at the dipole position by putting a unit dipole orientated along $j$ direction--can be obtain using numerical methods like boundary element method (BEM) or analytically.

    
What we are interested in is the relative ratio between the modified decay rates for a given $e\rightarrow g$ decay transition for all $ g $ and the natural linewidth, $ \Gamma_0 $, of the atoms, that is $ \Gamma_e^q/\Gamma_0 $. 
We can define an equivalent classical dipole, $ \mathbf{d}=\bra{g}\hat{\mathbf{d}}\ket{e} $, corresponding the $ e\rightarrow g $ decay transition which is along the $ \mathbf{e}_q $ ($q=\{\pm,0\}$) unit vector direction in the spherical irreducible harmonic basis, where
\begin{align}
    \mathbf{e}_\pm &=\mp \frac{\mathbf{e}_{\tilde{x}}\pm i\mathbf{e}_{\tilde{y}}}{\sqrt{2}}\\
    \mathbf{e}_0 &=\mathbf{e}_{\tilde{z}}
\end{align}
correspond to the $\sigma_\pm$ and $\pi$ transitions, respectively.
We have defined $ (\tilde{x},\tilde{y},\tilde{z}) $ as the quantization basis.
In our simulation of the square waveguide case, we have used $\varepsilon=4$ (without loss) to calculate the radiative mode induced Green's function tensor elements.
Therefore, the radiative mode caused decay rates can be calculated by
\begin{align}
\frac{\Gamma_{e,rad}^{q}}{\Gamma_0} &= 1+ \frac{\mathrm{Im}\left[\mathbf{e}_q\cdot \tensor{\mathbf{G}}_{ind,rad}(\mathbf{r}',\mathbf{r}')\cdot \mathbf{e}_q^*\right]}{ \mathrm{Im}\left[\mathbf{e}_q\cdot \tensor{\mathbf{G}}_0(\mathbf{r}',\mathbf{r}')\cdot \mathbf{e}_q^*\right]}
=1+ \frac{\mathbf{e}_q\cdot \mathrm{Im}\left[\tensor{\mathbf{G}}_{ind,rad}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{G_0(\mathbf{r}',\mathbf{r}')}.\label{eq:Gammaeq_rad}
\end{align}
Similarly, the guided mode contribution to the decay rates can be given by
\begin{align}
\frac{\Gamma_{e,gyd}^{q}}{\Gamma_0} &= \frac{\mathrm{Im}\left[\mathbf{e}_q\cdot \tensor{\mathbf{G}}_{gyd}(\mathbf{r}',\mathbf{r}')\cdot \mathbf{e}_q^*\right]}{ \mathrm{Im}\left[\mathbf{e}_q\cdot \tensor{\mathbf{G}}_0(\mathbf{r}',\mathbf{r}')\cdot \mathbf{e}_q^*\right]}
= \frac{\mathbf{e}_q\cdot \mathrm{Im}\left[ \tensor{\mathbf{G}}_{gyd}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{ G_0(\mathbf{r}',\mathbf{r}')}.\label{eq:Gammaeq_guide}
\end{align}
As has been derived in Ref~\cite{Qi2016},
\begin{align}
\mathrm{Im}[\tensor{\mathbf{G}}_{gyd}(\br',\br')] &= \pi \frac{\omega_{eg}}{v_g } \sum_{b, p} 
		\mathbf{u}_{b, p} (\br_{\!\perp}^\prime)\mathbf{u}^*_{b , p} (\br_{\!\perp}^\prime),
\end{align}
where $\mathbf{u}_{b, p} (\br_{\!\perp}^\prime)$ is the guided mode with propagation direction $ b=\pm 1 $ and polarization degeneracy index $ p $ at the dipole position in the transverse plane of the waveguide crossing section perpendicular to the waveguide axis.
$v_g$ is the group velocity of the degenerate guided modes.
For the nanofiber case, $ p=\pm 1 $ correspond to the right- and left-circularly polarized fundamental $\mathrm{HE}_{11}$ modes; 
for the \SWG case, $ p=\mathrm{H}/\mathrm{V} $ corresponding to the H- and V-modes chosen for the spin squeezing protocol.

\begin{figure}[htb]
\centering
 \begin{minipage}[h]{0.48\linewidth}
 %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
  \subfloat[h][$ \Gamma_{\rm nanofiber}/\Gamma_0 $]{
    \input{fig/nanofiber_decayrates.tex}
    %\includegraphics{fig/nanofiber_decayrates}
    \label{fig:nanofiber_decayrates}
    }
   \end{minipage}\vfill
   \begin{minipage}[h]{0.48\linewidth}
    %\begin{tabular}{*{2}{b{0.2\textwidth-2\tabcolsep}}}
     \subfloat[h][$\Gamma_{\rm SWG}/\Gamma_0$]{
       \input{fig/swg_decayrates.tex}
       %\includegraphics{fig/swg_decayrates}
       \label{fig:swg_decayrates}
       }
   \end{minipage}
   %\end{tabular}
\caption{Polarization-dependent decay rates for the nanofiber \protect\subref*{fig:nanofiber_decayrates} and SWG \protect\subref*{fig:swg_decayrates} as functions of the radial distance to the waveguides' symmetric centers. Line types differ the contributions from the guided and radiative modes and the combined total value; colors distinguish the polarizations of $ \sigma_\pm $, $ \pi $ and the averaged case. }\label{fig:decayrates}
\end{figure}

To calculate the modified decay rates of a multilevel atom in presence of the nanophotonic waveguides, one need to represent the dipole operators considering the particular level structure of the atom. 
In our case, we consider an alkali atom with the excited state $ \ket{e}=\ket{j'f'm'} $ indicated by the orbital angular momentum quantum number $ j' $, the total atomic angular momentum quantum number $ f' $ and the projected total atomic angular momentum quantum number $ m' $ of the hyperfine sublevel of the excited state; correspondingly, the ground state is indicated by $ \ket{g}=\ket{jfm} $. The dipole moment operator matrix elements defining the quantum transition from hyperfine sublevel $ m $ to $ m' $ can be given by
\begin{align}
\hat{\mathbf{d}}_{m'm} &= \sum_q \bra{j'f'm'}d_{f'f}^q\ket{jfm}\ket{j'f'm'}\bra{jfm}\mathbf{e}_q^*\\
&= \bra{f'}\vert d\vert\ket{f}\sum_q C_{f'm'}^{f,m;1,q}\ket{f'm'}\bra{fm}\mathbf{e}_q^*,
\end{align} 
where we have used the Wigner-Echart theorem to pull out the reduced matrix element $ \bra{f'}\vert d\vert\ket{f}=\bra{j'}\vert d\vert \ket{j}o_{jf}^{j'f'} $ with the relative oscillator strength $ o_{jf}^{j'f'} $ for a given nuclear angular moment quantum number; we have also denoted the Clebsch-Gordan coefficients $ C_{f',m'}^{f,m;1,q}=\bra{fm;1q}f'm'\rangle $~\cite{Deutsch2010a}.
Based on the selection rule, the Clebsch-Gordan coefficient $C_{f',m'}^{f,m;1,q}  $ is non-zero only if $ q+m=m' $.
Similarly, we can define the dipole moment operator defining the quantum transitions from hyperfine level $ f $ to $ f' $ by 
\begin{align}
\hat{\mathbf{d}}_{f'f} &= \sum_{m,m'} \hat{\mathbf{d}}_{m'm},
\end{align}
and from fine structure level $ j $ to $ j' $ by
\begin{align}
\hat{\mathbf{d}}_{j'j} &= \sum_{f,f'} \hat{\mathbf{d}}_{f'f}.
\end{align}
Given that the natural linewidth of the alkali atoms is determined by the quantum transitions between fine structure levels, which is proportional to the reduced dipole moment matrix element $ \bra{j'}\vert d\vert \ket{j} $, one can normalize the modified decay rates for particular levels to the natural linewidth by factorizing out the factor of $ \bra{j'}\vert d\vert \ket{j} $. For example, using the expressions of dipole moment operators, Eqs.~\eqref{eq:Gamma_e}, ~\eqref{eq:Gammaeq_rad} and~\eqref{eq:Gammaeq_guide} yield
\begin{align}
\frac{\Gamma_{m'\!,rad}}{\Gamma_0} &= 1+ \sum_f |o_{jf}^{j'f'}|^2\sum_{m,q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2\frac{\mathbf{e}_q\cdot \mathrm{Im}\left[\tensor{\mathbf{G}}_{ind,rad}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{G_0(\mathbf{r}',\mathbf{r}')}\\
\frac{\Gamma_{m'\!,gyd}}{\Gamma_0} &= \sum_f |o_{jf}^{j'f'}|^2\sum_{m,q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2 \frac{\mathbf{e}_q\cdot \mathrm{Im}\left[ \tensor{\mathbf{G}}_{gyd}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{ G_0(\mathbf{r}',\mathbf{r}')}
\end{align}
and 
\begin{align}
\frac{\Gamma_{f'\!,rad}}{\Gamma_0} &= 1+ \frac{1}{2f'+1}\sum_f |o_{jf}^{j'f'}|^2\sum_{m',m,q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2\frac{\mathbf{e}_q\cdot \mathrm{Im}\left[\tensor{\mathbf{G}}_{ind,rad}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{G_0(\mathbf{r}',\mathbf{r}')}\\
\frac{\Gamma_{f'\!,gyd}}{\Gamma_0} &= \frac{1}{2f'+1}\sum_f |o_{jf}^{j'f'}|^2\sum_{m',m,q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2 \frac{\mathbf{e}_q\cdot \mathrm{Im}\left[ \tensor{\mathbf{G}}_{gyd}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{ G_0(\mathbf{r}',\mathbf{r}')},
\end{align}
where the last two equations calculate the averaged decay rates for the $ f' $ sublevels due to the radiative and guided modes. 
In general, the dielectric modified atomic decay rate due to the transition $ \ket{f'm'} \rightarrow \ket{fm}$ can be given by 
\begin{align}
\frac{\Gamma_{m'\rightarrow m\!,rad}}{\Gamma_0} &= 1+ |o_{jf}^{j'f'}|^2\sum_{q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2\frac{\mathbf{e}_q\cdot \mathrm{Im}\left[\tensor{\mathbf{G}}_{ind,rad}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{G_0(\mathbf{r}',\mathbf{r}')}\\
\frac{\Gamma_{m'\rightarrow m\!,gyd}}{\Gamma_0} &= |o_{jf}^{j'f'}|^2\sum_{q=m'-m} \left|C_{f'm'}^{f,m;1,q}\right|^2 \frac{\mathbf{e}_q\cdot \mathrm{Im}\left[ \tensor{\mathbf{G}}_{gyd}(\mathbf{r}',\mathbf{r}')\right]\cdot \mathbf{e}_q^*}{ G_0(\mathbf{r}',\mathbf{r}')}.
\end{align}
For the decay rate equations above, whenever we sum over $ f $-levels, we should notice the frequency of the photon emission for different $ f $-levels should be different; given the hyperfine energy splitting ($ \sim 100 $MHz) and the wavelength of the light ($ \sim 1000 $nm), the frequency difference of the hyperfine level splitting is negligible. 

Due to the photon reflection from the dielectric interface, the frequency of the emitted light from the $ e\rightarrow g $ transition will also be shifted from the atomic resonance by 
\begin{align}
\delta\omega_{eg}=\omega-\omega_{eg}= \frac{1}{\hbar} \bra{g}\hat{\mathbf{d}}\ket{e}\cdot \mathrm{Re}\left[ \tensor{\mathbf{G}}_{ind,rad}(\br',\br';\omega_{eg})+\tensor{\mathbf{G}}_{gyd}(\br',\br';\omega_{eg})\right]\cdot \bra{e}\hat{\mathbf{d}}\ket{g}.
\end{align}
By noting the scattered electrical field $ \mathbf{E}_{scatt}(\br)=\tensor{\mathbf{G}}(\br,\br')\cdot \bra{e}\hat{\mathbf{d}}\ket{g} $, one can rewrite the normalized frequency shift $ \delta\omega_{m'm} $ to be
\begin{align}
\frac{\delta\omega_{m'm}}{\Gamma_0} &= \frac{1}{2} \frac{\left|o_{jf}^{j'f'}C_{f'm'}^{f,m;1,q=m'-m}\right|^2 \mathrm{Re}\left[ \mathbf{e}_q\cdot\left(\tensor{\mathbf{G}}_{ind,rad}(\br')+\tensor{\mathbf{G}}_{ind,gyd}(\br')\right)\cdot\mathbf{e}_q^*\right]}{G_0(\br',\br')}\\
&=\frac{1}{2} \frac{\left|o_{jf}^{j'f'}C_{f'm'}^{f,m;1,q=m'-m}\right|^2 \mathrm{Re}\left[ \mathbf{e}_q\cdot\left(\mathbf{E}_{ind,rad}^{(q^*)}(\br')+\mathbf{E}_{ind,gyd}^{(q^*)}(\br')\right)\right]}{G_0(\br',\br')},
\end{align}
where $ \mathbf{E}_{ind,rad/gyd}^{(q^*)}(\br') $ is the induced radiative/guided mode component of the local electrical field due to a dipole polarized along $ \mathbf{e}_q^* $ direction, which can be calculated using the BEM simulation. 
By definition, this is the Lamb shift due to the light scattering from the dielectric interface and has excluded the propagating field in vacuum. 

In another way, the Lamb shift can also be calculated through evaluating the real part of the medium's Green's function tensors from their imaginary parts using the Kramers-Kronig relations~\cite{Dzsotjan2011}. In that case, the total Green's function tensor and the vacuum Green's function tensor may be required to be evaluated accordingly, but the approximate solution for a cylindrical waveguide geometry has been given in Ref.~\cite{Dzsotjan2011} which works pretty well when the atoms are not too close to the waveguide surface as the Lamb shift decays on the order of $ 1/(r_\perp-r_\perp^0)^3 $ where $ r_\perp^0 $ is the radial coordinate of the waveguide surface. One can use the approximate solution for a cylindrical waveguide to validate the accuracy of the Lamb shift calculation using the BEM approach. 

 
The collective spin dynamics in the QND measurement and spin squeezing process is described by the stochastic master equation defined in Eq.~\eqref{eq:totaldrhodt}. The optical pumping dynamics of the $ j $-th atom are governed by 
\begin{align}
\left.\dt{\hat{\rho}^{(j)}}\right|_{op} &= \mathcal{D}[\hat{\rho}^{(j)}]=-\frac{i}{\hbar}\left\{\hat{H}_{\rm eff}^{(j)},\hat{\rho}^{(j)} \right\}_+ + \gamma_s\sum_{q}\hat{W}_q(\br'_j)\hat{\rho}^{(j)}\hat{W}_q^\dagger(\br'_j) %\\
%&=-\gamma_s\frac{1}{2}\sum_{a,b}\gamma_{ba}\ket{b}\bra{a}+
%\!\!\!\!\!\!\sum_{q,q',q'',a,b,c,d,f',f'',m',m''}\!\!\!\!\!\! \gamma_sw_{dcq}^{f''m''q''}\left(w_{abq}^{f'm'q'}\right)^*\ket{d}\bra{c}\hat{\rho}^{(i)}\ket{b}\bra{a}.
\end{align} 
We have defined a characteristic photon scattering rate, $\gamma_s \equiv \frac{\Gamma_0\Omega^2}{4\Delta_{F}^2}= \frac{\sigma_0}{A_{\rm in}}\frac{\Gamma_0^2}{4 \Delta_{F}^2} \dot{N}_L $ with an effective detuning $ \Delta_F $ defined by $ \frac{1}{\Delta_F}=\sum_{f'}\frac{C_{ff'}^{(1)}}{\Delta_{ff'}} $ and $ \Delta_{ff'}=\omega-\omega_{ff'} $.
We have also defined Rabi frequency $ \Omega=2\bra{j}|d|\ket{j'}\mathcal{E}^{(+)}_{\rm in}/\hbar $ with reduced optical dipole matrix element $\bra{j}|d|\ket{j'}$ and field amplitude $ \mathcal{E}^{(+)}_{\rm in}=|\mathbf{E}_{\rm in}^{(+)}(\br')| $.
The total rate of photon scattering by an atom from the $\ket{a}\equiv \ket{f,f_x=a}$ to $ \ket{b}\equiv\ket{f,f_x=b} $ hyperfine ground state in the $x$-basis is
	\begin{equation}\label{Eq::gammaf}
		\gamma_{ba}=- \frac{2}{\hbar} {\rm Im} \big[ \bra{f,b} \hat{h}_{\rm eff}\ket{f,a} \big] ,
	\end{equation}
The effective nonHermitian light-shift Hamiltonian for one atom (labeled with superscript $ (i) $) is given by
\begin{align}
\hat{H}_{\rm eff}^{(i)} \equiv \hat{h}_{\rm eff}&= - \hat{\mathbf{E}}^{(-)}_{\rm in}(\mathbf{r}' ; t ) \cdot \poltens \cdot \hat{\mathbf{E}}^{(+)}_{\rm in}(\mathbf{r}' ;t ),
\end{align}
where the polarizability operator $  \poltens=\sum_{f',q}\hat{\tensor{\mathbf{A}}}(f,f',q)$ with elements of $ \hat{\tensor{\mathbf{A}}} $ given by
\begin{align} \label{Eq::PolarizabilityIrrep}
		\hat{A}_{ij}(f,f',q)&\equiv -\frac{\sigma_0}{8\pi k_0\gamma_s}\frac{\Gamma_{f'\!\!,\, i}^q}{\Delta_{ff'}^q+i\Gamma_{f'\!\!,\, i}^q/2}\hat{e}_i^*\cdot\hat{\mathbf{D}}_{ff'}\hat{\mathbf{D}}_{f'f}^\dagger \cdot \hat{e}_j \\
		&=  -\frac{\sigma_0}{8\pi k_0\gamma_s}\frac{\Gamma_{f'\!\!,\, i}^q}{\Delta_{ff'}^q\!+\! i\Gamma_{f'\!\!,\, i}^q/2}\left\{ C_{ff'}^{(0)} \delta_{i,j}\hat{\mathbbm{1}}\!+\! iC_{ff'}^{(1)}\epsilon_{ijk}\hat{f}_k \!+\! C_{ff'}^{(2)} \Big[ \smallfrac{1}{2} ( \hat{f}_i\hat{f}_j \!+\!\hat{f}_j\hat{f}_i )\!-\!\smallfrac{1}{3} \hat{\mathbf{f}}\!\cdot\!\hat{\mathbf{f}} \delta_{i,j} \Big]\right\}, 
\end{align}
where $\hat{\mathbf{f}}$ is the atomic spin operator in hyperfine multiplet $f$, and $ \epsilon_{ijk} $ is the Levi-Civita symbol. 
Definitions of the interaction coefficients $ C_{ff'}^{(n)} $ can be found in Ref.~\cite{Qi2016}, which correspond to scalar, vector and tensor atom-light interactions for $ n=0,1,2 $, respectively.

Given the geometry of the Faraday spin squeezing protocol for optical nanofiber and square waveguides discussed in this paper, the local electric field at the atom positions is linearly polarized. 
By denoting the local field's polarization direction as the $ x$-direction and the propagation direction of the guided light as the $ z $-direction, the effective atom-light interaction Hamiltonian in a static reference frame can be written as 
\begin{align}
\hat{h}_{\rm eff} &= -\frac{i\hbar}{2}\sum_{f'} \gamma'_s \left[C_{ff'}^{(0)}\hat{\mathbbm{1}}+C_{ff'}^{(2)}(\hat{f}_{x}^2-\frac{\hat{\mathbf{f}}^2}{3} ) \right],\label{eq:Heff_Faraday_C02}
\end{align}
where the intrinsic photon scattering rate $ \gamma'_s\equiv \frac{\Gamma_{f'}^x\Omega^2}{4(\Delta_{ff'}^2+\left(\Gamma_{f'}^x\right)^2/4 )}$ and in the far-detuning regime, $\gamma'_s \approx \frac{\Gamma_{f'}^x\Omega^2}{4\Delta_{\rm eff}^2}=\frac{\sigma_0}{\Ain}\left(\frac{\Gamma_{f'}^x}{2\Delta_{\rm eff}} \right)^2\dot{N}_L $. Note here the vector interaction term vanishes given a linearly polarized light; we have ignored the energy shift, which is valid when the detuning is much larger than the hyperfine level splitting.
In our simulations of spin squeezing in this paper, we have used the decay rate $ \Gamma_{f',q} $ to indicate the decay rates from the hyperfine $ f' $ manifold sublevels of excited states with a photon emission polarized along the $ \mathbf{e}_q $ direction; the effective detuning is also an averaged detuning from the fine structure excited level $ j' $ to the ground fine structure manifold $ j $ with resonant frequency $ \omega_D $ of $ D_1 $ line transitions--that is $ \Delta_{\rm eff}=\omega -\omega_D $ with probe frequency at $ \omega $ in vacuum given a far-detuned $ \sim 1 $nm of detuning from the $ D_1 $ line transition of $ ^{133}Cs $ atoms. 
Compared to the normal characteristic photon scattering rate $ \gamma_s $, we can see they are defined in different scales.
In general, they are related given a transition between ground hyperfine structure level $ f $ and excited hyperfine structure level $ f' $ by 
\begin{align}
\gamma'_s(f')=\gamma_s \frac{\Delta_F^2}{\Delta_{ff'}^2},
\end{align}
and hence $ \frac{\sqrt{\Gamma_{f'}^x}\Omega/2}{\Delta_{ff'}\pm i\Gamma_{f'}^x/2}\approx \frac{\sqrt{\Gamma_{f'}^x}\Omega/2}{\Delta_{ff'}}=\sqrt{\gamma_s}\frac{\Delta_F}{\Delta_{ff'}} $ in the far-detuning regime.

We define the Lindblad jump operators of optical pumping among ground states by~\cite{Deutsch2010a}
	\begin{align}\label{Eq::Wq_Faraday}
		\hat{W}_q &= \frac{1}{\sqrt{\gamma_s}}\sum_{f'}\frac{\sqrt{\Gamma_{f'}^q}\Omega/2}{\Delta_{f'f}^q+i\Gamma_{f'}^q/2}\mathbf{e}_q^*\cdot(\hat{\mathbf{D}}_{ff'}  \hat{\mathbf{D}}^\dagger_{f'f} )\cdot\mathbf{e}_{\rm in} \\
		&= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'k}\! \frac{\sqrt{\Gamma_{f'}^q}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^q/2} \left[\delta_{qx}C_{ff'}^{(0)}\hat{\mathbbm{1}} \!+\! iC_{ff'}^{(1)}\epsilon_{qxk}\hat{f}_k \!+\! C_{ff'}^{(2)}\left(\frac{\hat{f}_q\!\hat{f}_{x}\!+\!\hat{f}_{x}\!\hat{f}_q }{2} \!-\! \frac{\delta_{qx}}{3}\hat{\mathbf{f}}^2 \right) \right].
%		&=\sum_{f',m',q',a,b} w_{baq }^{f'm'q'}\ket{b}\bra{a}.
	\end{align}
Each jump operator $\hat{W}_q$ is associated with absorption of the probe photon polarized along $ \mathbf{e}_{\rm in} $ followed by spontaneous emission of a photon with polarization $ \mathbf{e}_q $, where $q= \{0,\pm 1\}$ labels spherical basis elements for $\pi$ and $ \sigma_\pm$ transitions. 
%Here the dimensionless raising operator $ \mathbf{e}_q\cdot\hat{\mathbf{D}}_{f'f}^\dagger= \sum_{m',m} o_{jf}^{j'f'} C_{f',m'}^{f,m;1, q}\ket{f',m'}\bra{f,m} $,
%where $ C_{f',m'}^{f,m;1, q}=0 $ unless $ m'=m+q $ with $ C_{f',m+q}^{f,m;1, q}=\Braket{f',m+q}{f,m;1,q}$ being the Clebsch-Gordan coefficients, and
%\begin{equation}
%\big| o_{jf}^{j'f'} \big|^2=(2j'+1)(2f+2) \bigg\{
%\begin{array}{ccc}
%f' & 7/2 & j' \\
% j & 1 & f
% \end{array}
% \bigg\}
%\end{equation}
%are the relative oscillator strengths determined by the relevant Wigner 6-$J$ symbol.
%In our protocols, we assume the probe light is so far-detuned from any of the atomic resonances that the tilting of the hyperfine structure levels due to an external magnetic field becomes irrelevant and we can set $ \Delta_{ff'}^q=\Delta_{ff'} $ and $ \Delta_{ff'}\gg \Gamma_f'^q $ for arbitrary $ f' $ and $ q $.

Therefore, in the static $ \left\{x,y,z \right\} $ basis, we have
\begin{align}
\hat{W}_{x} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^x}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^x/2} \left[C_{ff'}^{(0)}\hat{\mathbbm{1}} + C_{ff'}^{(2)}\left(\hat{f}_{x}^2-\frac{1}{3}\hat{\mathbf{f}}^2 \right) \right]\\
\hat{W}_{y} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^y}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^y/2} \left(-iC_{ff'}^{(1)}\hat{f}_z + C_{ff'}^{(2)}\frac{\hat{f}_{x}\hat{f}_{y}+\hat{f}_{y}\hat{f}_{x}}{2} \right)\\
\hat{W}_{z} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^z}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^z/2} \left(iC_{ff'}^{(1)}\hat{f}_{y} + C_{ff'}^{(2)}\frac{\hat{f}_z\hat{f}_{x}+\hat{f}_{x}\hat{f}_z}{2}  \right).
\end{align}


The optical pumping dynamics of the $ j $-th atom are governed by 
\begin{align}
\left.\dt{\hat{\rho}^{(j)}}\right|_{op} &= \gamma_s\mathcal{D}[\hat{\rho}^{(j)}]=-\frac{i\gamma_s}{\hbar}\left\{\hat{h}^{\rm eff}_{\rm eff},\hat{\rho}^{(j)} \right\}_+ + \gamma_s\sum_{q}\hat{W}_q(\br'_j)\hat{\rho}^{(j)}\hat{W}_q^\dagger(\br'_j)\\
\\
&=-\gamma_s\frac{1}{2}\sum_{a,b}\gamma_{ba}\ket{b}\bra{a}+
\!\!\!\!\!\!\sum_{q,q',q'',a,b,c,d,f',f'',m',m''}\!\!\!\!\!\! \gamma_sw_{dcq}^{f''m''q''}\left(w_{abq}^{f'm'q'}\right)^*\ket{d}\bra{c}\hat{\rho}^{(i)}\ket{b}\bra{a}
\end{align}
Eqs.~\eqref{Eq::gammaf} and~\eqref{Eq::Wq_Faraday} yield,
\begin{subequations}
	\begin{align}
		\gamma_{ba} 
		&=\frac{n_g\dot{N}_L}{\gamma_s}  \sum_{f',q} \sigma (\Delta_{ff'} ) \mathbf{u}^*_\inp(\br'_\perp)\cdot \bra{b} \hat{\tensor{\mbf{A}}}(f,f') \ket{a}  \cdot \mathbf{u}_\inp(\br'_\perp)\\
		&\approx  \sum_{f',m'} \frac{\Delta_{F}^2}{\Delta_{ff'}^2}\sum_{q,q'} \big| o_{jf}^{j'f'} \big|^2C_{f',b+q'}^{f,b;1, q'}C_{f',a+q}^{f,a;1, q} \mathbf{e}_{q'}^* \cdot (\mathbf{e}_{\rm in}\mathbf{e}_{\rm in}^* )\cdot \mathbf{e}_q,
	\end{align}
\end{subequations}
	\begin{align}
		w_{baq}^{f'm'q'}
		&\approx  \frac{\Delta_{F}}{\Delta_{ff'}+i\Gamma_{f'}^q/2} \big| o_{jf}^{j'f'}  \big|^2 C_{f'm'}^{f,b;1 q}C_{f',m'}^{f,a;1,q'} (\mathbf{e}_{q'}^* \cdot \mathbf{e}_{\rm in}),
	\end{align}
where $ \sigma (\Delta_{ff'} )  = \sigma_0 \Gamma_0^2/4\Delta^2_{f' f}$ is the the scattering cross section at the probe detuning in free space. 

Now, we consider a static magnetic field is applied to the atoms to fix the quantization axis to the $ \mathbf{e}_{z} $ direction pointing along the waveguide axis.
We assume the magnetic field is so strong that the Larmor processing is much faster than the atomic decay and atom-photon interaction processes and the transverse components of the atomic angular momentum operators will be averaged out in the process of spin squeezing dynamics.
In theory, this leads us to transfer the master equations of the collective spin dynamics to the rotating frame determined by the fast-rotating transform operator
\begin{align}
\hat{U}_B(t) &= e^{-i\Omega_Bt\hat{f}_z},
\end{align}
where $ \Omega_B $ is the Larmor processing frequency of the external magnetic field.
In the rotating frame, a quantum operator $ \hat{A} $ is transfered into $ \hat{A}' $ through $ \hat{A}\rightarrow \hat{A}'=\expect{\hat{U}_B^\dagger\hat{A}\hat{U}_B }_T $, where the notation $ \expect{\cdot}_T=\frac{1}{T}\int\cdot dt $ is the time average of observables in a period $ T $.
The density operator preserves its form in the rotating frame.
We can solve the transformed master equations by employing the Baker-Campbell-Hausdorff formula that $ e^{\lambda\hat{A}}\hat{B}e^{-\lambda\hat{A}}=\sum_{n=0}^\infty\frac{\lambda^n}{n!}\hat{C}_n $, where $ \hat{C}_0=\hat{B} $ and $ \hat{C}_n=\left[\hat{A},\hat{C}_{n-1} \right] $ for $ n>1 $, and the commutators of atomic angular momentum operators, $ \left[\hat{f}_m, \hat{f}_n\right]=i\sum_p\epsilon_{mnp}\hat{f}_p $.
The following static-rotating frame transformation relationships can be proved easily:
\begin{subequations}\label{eq:rotationtransf}
	\begin{align}
	\hat{U}_B^\dagger\hat{f}_{x}\hat{U}_B&=\cos(\Omega_Bt)\hat{f}_x-\sin(\Omega_Bt)\hat{f}_y,\\
	\hat{U}_B^\dagger\hat{f}_{y}\hat{U}_B &= \sin(\Omega_Bt)\hat{f}_x+ \cos(\Omega_Bt)\hat{f}_y, \\ \hat{U}_B^\dagger\hat{f}_{z}\hat{U}_B &=\hat{f}_z,\quad \hat{U}_B^\dagger\hat{\mathbbm{1}}\hat{U}_B =\hat{\mathbbm{1}}.
	\end{align}
\end{subequations}
In the rotating frame, operators with a transverse atomic angular momentum will be averaged to vanish. For example,
\begin{subequations}\label{eq:rotationtransf_f}
	\begin{align}
	\hat{f}_{x}&=\hat{f}_{y} \rightarrow 0, \quad \hat{f}_{z}\rightarrow\hat{f}_z, \quad \hat{f}_{z}^ 2\rightarrow\hat{f}_z^2,\\
	\hat{f}^2_{x} &= \hat{f}^ 2_{y} \rightarrow \frac{1}{2}(\hat{\mathbf{f}}^2-\hat{f}_z^2),\\
	\hat{f}_{x}\hat{f}_{y} &\rightarrow\frac{1}{2}\hat{f}_z,\quad \hat{f}_{y}\hat{f}_{x}\rightarrow -\frac{1}{2}\hat{f}_z,\quad \hat{f}_{i=x,y}\hat{f}_{z}\rightarrow 0.
	\end{align}
\end{subequations}



Using the transformation relationships defined in Eqs.\eqref{eq:rotationtransf}, the loss Hamiltonian in the rotating frame becomes
\begin{subequations}\label{eq:rotationtransf_hloss}
\begin{align}
\hat{h}_{\rm eff} =-\frac{i\hbar}{2} \sum_{f'}\gamma'_s \left[C_{ff'}^{(0)}\hat{\mathbbm{1}} + \frac{C_{ff'}^{(2)}}{6}(\hat{\mathbf{f}}^2-3\hat{f}_z^2 ) \right],
\end{align}
\end{subequations}
where $ z $-direction is the waveguide axis direction, and $ \hat{\mathbf{f}}^2=\hat{\mathbf{f}}\cdot\hat{\mathbf{f}} $.


By using the transformation relationships of Eqs.\eqref{eq:rotationtransf}, the jump operators become 
\begin{subequations}\label{eq:rotationtransf_Wxyz}
\begin{align}
\hat{W}_{x} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^x}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}/2} \left[C_{ff'}^{(0)}\hat{\mathbbm{1}} + \frac{C_{ff'}^{(2)}}{6}\left(\hat{\mathbf{f}}^2-3\hat{f}_{z}^2 \right) \right]\\
\hat{W}_{y} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} -\frac{i\sqrt{\Gamma_{f'}^y}\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^y/2} C_{ff'}^{(1)}\hat{f}_z  \\
\hat{W}_{z} &= 0.
\end{align}
\end{subequations}

By using the fact that 
\begin{subequations}
\begin{align}
\hat{\mathbf{f}}^2 &=f(f+1)\hat{\mathbbm{1}}\\
\hat{f}_z &=\sum_{m=1}^{2f+1}(f-m+1)\hat{\sigma}_{mm}\\
\hat{f}_z^2 &= \sum_{m=1}^{2f+1}(f-m+1)^2\hat{\sigma}_{mm}
\end{align}
\end{subequations}
in the rotating frame, both $ \hat{h}_{\rm eff} $ and $ \hat{W}_{q} $ become diagonal, and Eqs.\eqref{eq:rotationtransf_hloss} and~\eqref{eq:rotationtransf_Wxyz} can be simplified as
\begin{subequations}
\begin{align}
\hat{h}_{\rm eff} &= -\frac{i\hbar}{2} \sum_{f'} \gamma'_s(f') \sum_{m=-f}^f\left[C_{ff'}^{(0)} + \frac{C_{ff'}^{(2) }}{6}(f(f+1)-3m^2) \right]\hat{\sigma}_{mm}\\
\hat{W}_{x} &= \frac{1}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^x }\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^x/2 } \sum_{m=-f}^f\left[C_{ff'}^{(0)} + \frac{C_{ff'}^{(2) }}{6}(f(f+1)-3m^2) \right]\hat{\sigma}_{mm}\\
\hat{W}_{y} &= -\frac{i}{\sqrt{\gamma_s}}\!\sum_{f'} \frac{\sqrt{\Gamma_{f'}^y }\Omega/2}{\Delta_{ff'}+i\Gamma_{f'}^y/2 } \sum_{m=-f}^f C_{f'ff'}^{(1)}m\hat{\sigma}_{mm}\\
\hat{W}_z &=0.
\end{align}
\end{subequations}
After some algebra, one can obtain the optical pumping master equation in the rotating frame as
\begin{align}
\left. \dt{\hat{\rho}}\right|_{\rm op} 
&= - \sum_{f'} \gamma'_s(f') \left[\left(C_{ff'}^{(0)}+\frac{f(f+1)}{12}C_{ff'}^{(2)} \right)\hat{\rho}-\frac{C_{ff'}^{(2)}}{4}(\hat{f}_z^2\hat{\rho}+\hat{\rho}\hat{f}_z^2) \right]\nn\\
&\quad+\sum_{f',f''} \frac{\sqrt{\Gamma_{f'}\Gamma_{f''} }\Omega^2/4 }{\Delta_{ff'}\Delta_{ff''}+\Gamma_{f'}\Gamma_{f''}/4+i(\Delta_{ff''}\Gamma_{f'}-\Delta_{ff'}\Gamma_{f''} ) }\nn\\
&\quad\cdot\left\{C_{ff'}^{(0)}C_{ff''}^{(0)}\hat{\rho}+ C_{ff'}^{(0)}C_{ff''}^{(2)}\frac{\hat{\rho}}{6}(\fo^2-3\hat{f}_z^2) + C_{ff''}^{(0)}C_{ff'}^{(2)}(\fo^2-3\hat{f}_z^2)\frac{\hat{\rho}}{6} \right.\nn\\
&\quad\quad + \frac{1}{2}C_{ff'}^{(1)}C_{ff''}^{(1)}(\hat{f}_x\hat{\rho}\hat{f}_x+\hat{f}_y\hat{\rho}\hat{f}_y+2\hat{f}_z\hat{\rho}\hat{f}_z )\nn\\
&\quad\quad -\frac{1}{4}C_{ff'}^{(1)}C_{ff''}^{(2)}(\fx\rhoo\fx-\fy\rhoo\fy-2i\fx\rhoo\fz\fy+2i\fy\rhoo\fx\fz )\nn\\
&\quad\quad +\frac{1}{4}C_{ff''}^{(1)}C_{ff'}^{(2)}(\fx\rhoo\fx-\fy\rhoo\fy-2i\fz\fy\rhoo\fx+2i\fx\fz\rhoo\fy ) \nn\\
&\quad +C_{ff'}^{(2)}C_{ff''}^{(2)}\left[\frac{1}{4}f^2(f+1)^2\rhoo-\frac{f(f+1)}{6}(\fo^2-\fz^2)\rhoo-\frac{f(f+1)}{6}\rhoo(\fo^2-\fz^2) \right.\nn\\
&\quad\quad\quad+\frac{1}{2}\fx^2\rhoo\fx^2+\frac{1}{2}\fy^2\rhoo\fy^2+\frac{1}{4}(i\fz+2\fy\fx)\rhoo(i\fz+2\fy\fx)\nn\\
&\left.\left.\quad\quad\quad +\frac{1}{8}(i\fy+2\fx\fz)\rhoo(i\fy+2\fx\fz)+\frac{1}{8}(i\fx+2\fz\fy)\rhoo(i\fx+2\fz\fy) \right]\right\}.
\end{align}
As a sanity check, we consider a far-detuning regime where the detuning on hyperfine sublevels is indistinguishable so that we can denote $ f''=f' $ and ignore all tensor polarizability terms where $ C_{ff'}^{(2)} $ presents, and the optical pumping part of the master equation above becomes
\begin{align}
\left.\dt{\rhoo}\right|_{\rm op} &= \gamma'_s  \left[C_{fj'}^{(0)}(C_{fj'}^{(0)}-1)\rhoo+\frac{1}{2}(C_{fj'}^{(1)})^2(\fx\rhoo\fx+\fy\rhoo\fy+2\fz\rhoo\fz ) \right]\\
&= -\frac{2\gamma'_s}{9}+\frac{\gamma'_s}{18f^2}(\fx\rhoo\fx+\fy\rhoo\fy+2\fz\rhoo\fz )
\end{align}
for the far-detuned $ D_1 $- and $ D_2 $-line transitions of cesium atoms,
which is a well-known result in previous studies~\cite{Deutsch2010a,Baragiola2014}.
Above, we have defined $ C_{fj'}^{(0)}=\sum_{f'}C^{(0)}_{ff'}=1/3  $ for $ j'=1/2 $ or $ 2/3 $ for $ j'=3/2 $; $ C^{(1)}_{fj'}=\sum_{f'}C^{(1)}_{ff'}=\pm g_f/3 $ for $ j'=1/2 $ and $ j'=3/2 $, respectively, and $ g_f=1/f $ for our case. 


%===================APPENDIX: Equations of motion =====================%
\section{Generalized master equation for spin squeezing dynamics with arbitrary spin number $f$} \label{Appendix::OpticalPumpingForGeneralF}
\comment{Todo: see notes on QND measurement and spin squeezing using Faraday interactions (part VI).}


%===================APPENDIX: Relationships between collective spin operators and microscopic operators =====================%
\section{Relationships between some collective operators and the first- and second-moments in the symmetric subspace} \label{Appendix::collectivespinoperators}

\comment{Todo: Relations between collective angular momentum operators and microscopic $\sigma_{ba}$ operators in the qubit and qutrit symmetric subspace. See handwritten notes: QND measurement and spin squeezing using Faraday interactions (part IV).}
\begin{subequations}
	\begin{align}
	\expect{\hat{F}_x} &= \sum_i^{N_A}\expect{\hat{f}_x}\nonumber\\
	&= -N_A \left[f\expect{\sigmauu}+(f-1)\expect{\sigmadd}+(f-2)\expect{\sigmatt } \right]\label{eq:Fx_qutrit}\\
	\expect{\Delta F_z^2} &= N_A\expect{\Delta f_z^2} + N_A(N_A-1)\expect{\Delta f_z^{(1)}\Delta f_z^{(2)} }_s\nn\\
	&=N_A\left\{ \frac{f}{2}\left(\expect{\sigmauu}+\expect{\sigmadd}-2\expect{\sigmaud}\expect{\sigmadu}-\expect{\sigmaud}^2-\expect{\sigmadu}^2 \right)\right. \label{eq:DeltaFz2_fz}\\
	&\quad\quad+ \frac{2f-1}{2}\left(\expect{\sigmadd}+\expect{\sigmatt}-2\expect{\sigmadt}\expect{\sigmatd}-\expect{\sigmadt}^2-\expect{\sigmatd}^2 \right)\nn\\
	&\quad\quad +\left. \frac{\sqrt{f(2f-1)}}{2}\left[\expect{\sigmaut }+\expect{\sigmatu} -2\expect{\sigmaud}(\expect{\sigmadt}+\expect{\sigmatd}) -2\expect{\sigmadu}(\expect{\sigmadt}+\expect{\sigmatd} ) \right] \right\}\nn\\
	&\quad +N_A(N_A-1)\left[\frac{f}{2}(\expect{\Dsigmaud^{(1)}\Dsigmaud^{(2)} }_s +2\expect{\Dsigmaud^{(1)}\Dsigmadu^{(2)} }_s+\expect{\Dsigmadu^{(1)}\Dsigmadu^{(2)} }_s )\right.\label{eq:DeltaFz2_qutrit}\\
	&\quad\quad + \frac{2f-1}{2}(\expect{\Dsigmadt^{(1)}\Dsigmadt^{(2)} }_s +2\expect{\Dsigmadt^{(1)}\Dsigmatd^{(2)} }_s +\expect{\Dsigmatd^{(1)}\Dsigmatd^{(2)} }_s) \nn\\
	&\quad\quad + \left. \sqrt{f(2f-1)}(\expect{\Dsigmaud^{(1)}\Dsigmadt^{(2)} }_s +\expect{\Dsigmaud^{(1)}\Dsigmatd^{(2)} }_s +\expect{\Dsigmadu^{(1)}\Dsigmadt^{(2)} }_s+\expect{\Dsigmadu^{(1)}\Dsigmatd^{(2)} }_s )\right]\nn
	\end{align}
\end{subequations}

Certainly, the collective spin dynamics can also be calculated in the $ \left\{\left.\expect{\hat{f}_i},\expect{\Delta f_i^{(1)}\Delta f_j^{(2)} }_s\right|i,j=0,x,y,z \right\} $ basis with $ \hat{f}_0=\hat{\mathbbm{1}} $. For example, one can prove that, by ignoring couplings between $ \ket{\uparrow} $ and $ \ket{T} $, Eqs.~\ref{eq:Fx_qutrit} and~\ref{eq:DeltaFz2_qutrit} lead to 
\begin{subequations}
\begin{align}
\expect{\hat{F}_x} &= -f \expect{\hat{N}_\uparrow}-(f-1)\expect{\hat{N}_\downarrow}-(f-2)\expect{\hat{N}_T}\\
\expect{\Delta F_z^2} &\approx \expect{\hat{N}_f}\expect{\Delta f_z^2}+\expect{\hat{N}_f}(\expect{\hat{N}_f}-1)\expect{\Delta f_z^{(1)}\Delta f_z^{(2)} }\!_s,
\end{align}
\end{subequations}
where $ \hat{N}_i=\sum_n^{N_A}\hat{\sigma}_{ii}^{(n)} $ (for $ i=\uparrow,\downarrow,T $) and $ \hat{N}_f=\hat{N}_\uparrow+\hat{N}_\downarrow+\hat{N}_T $ are the population operators. Then the problem is to establish and solve the stochastic equations for the expectation values and covariances of the population and angular momentum operators. For example, QND measurement process which generates the spin squeezing dynamics is governed by
\begin{align}
d\expect{\Delta f_z^2}|_{QN\!D} = -\kappa \left[\expect{\Delta f_z^2}+(N_A-1)\expect{\Delta f_z^{(1)}\Delta f_z^{(2)} }_s \right]^2dt=d\expect{\Delta f_z^{(1)}\Delta f_z^{(2)} }_s|_{QN\!D}.
\end{align} 
It seems the spin squeezing dynamics is purely determined by the $ \expect{\Delta f_z^2} $ and $ \expect{\Delta f_z^{(1)}\Delta f_z^{(2)} }$, but in deriving the equations for the full set of operator basis considering both the optical pumping and QND measurement contributions, all operators become coupled to each other.


\end{appendix}

\end{document}
